Load libraries

Load data Name the data you export from FileMaker Pro by their exact table names and save them as CSVs, e.g. Larvae_Morphology.csv

#setwd("~/Repos/LarvaeTransGen2018")
setwd("~/R/GitHub/LarvaeTransGen2018/data") #Elise's working directory
#upload all of the data tables for the larvae experiment
Larvae_Morphology <- read.csv("../data/Larvae_Morphology.csv") #contains data from CellProfiler from the larvae outlines
Barcode_Jar <- read.csv("../data/Barcode_Jar.csv") #contains data on CrossID, seatable, and whether or not larvae were present
Block_ID <- read.csv("../data/Block_ID.csv") #contains data on the block
Cross<- read.csv("../data/Cross.csv") #contains data on the female and male IDs for the cross and whether the cross was for QG or Meth
Fert_QG<- read.csv("../data/Fert_QG.csv") #contains data on fertilization counts, the JarIDs for the crosses
Header_WC<- read.csv("../data/Header_WC.csv") #contains data on the header tanks used to fill the jars
Larvae_Calcification <- read.csv("../data/Larvae_Calcification.csv") #contains data on filter weights and counts
Larvae_Counts <- read.csv("../data/Larvae_Counts.csv") #contains data on sedgewick rafter larvae counts
Larvae_WC <- read.csv("../data/Larvae_WC.csv") #contains data on larvae jar water chemistry
Parent_ID <- read.csv("../data/Parent_ID.csv") #contains data on the adults at the time of shucking. Includes Parent ID assignments
Larvae_Cilia<- read.csv("../data/Larvae_Cilia.csv") #contains data on cilia extrusion of the larvae photographed for morphology
Egg_Morphology<- read.csv("../data/Egg_Morphology.csv") #contains data on egg sizes from CellProfiler
Adult_Sample<- read.csv("../data/Adult_Sample.csv") #contains data on adult oysters at the time of collection
WC_Standard<- read.csv("../data/WC_Standard.csv") #contains data on the standards used for water chemistry
Tank_WC<- read.csv("../data/Tank_WC.csv") #contains data on the adult tank water chemistry

Make fields factors:

Adult_Sample$AccTankID<- as.factor(Adult_Sample$AccTankID)
Adult_Sample$TrtTankID<- as.factor(Adult_Sample$TrtTankID)
Barcode_Jar$JarSeatable<- as.factor(Barcode_Jar$JarSeatable)
Barcode_Jar<- subset(Barcode_Jar, JarSeatable =="2" | JarSeatable =="4" | JarSeatable =="6")
Tank_WC$ParentTrt<- as.factor(Tank_WC$ParentTrt)
Tank_WC$Tank<- as.factor(Tank_WC$Tank)
Tank_WC$Folder<- as.factor(Tank_WC$Folder)
WC_Standard$CRM<- as.factor(WC_Standard$CRM)
Larvae_WC$DateFilter<- as.factor(Larvae_WC$DateFilter)

Calculate average larval jar water chemistry based on the subset of bottles that were run on the VINDTA.

#start by joining Barcode_Jar and Larvae_WC datasets 
JarChem<- merge(x=Barcode_Jar, y= Larvae_WC, by="JarID", all.x=TRUE)
#combine Date filter and JarTrt columns to get a unique combo value for each
JarChem<-unite(JarChem,BlockTrt, c("DateFilter","JarTrt"), sep='', remove=F)
#get only jars with larvae
JarChem<- subset(JarChem, Larvae =="TRUE")
JarChem<- subset(JarChem, DateFilter =="20180713")
boxplot(JarpHSW~BlockTrt, data=JarChem, ylim=c(6,8))

#get means for the parameters
JarMeans<- ddply(JarChem, .(BlockTrt), numcolwise(mean, na.rm=T))
#remove the unnecessary columns from JarMeans
JarMeans<- subset(JarMeans, select= c(BlockTrt,JarAlkCalc, JarpCO2, JarHCO3, JarCO3, JarCO2, JarOmegaCalcite, JarOmegaArg))   
#add Est for "estimated" to the column names
colnames(JarMeans)<- paste("Est", colnames(JarMeans), sep="")
#add columns to JarChem for mean carb values
JarChem<- merge(x=JarChem, y=JarMeans, by.x="BlockTrt", by.y="EstBlockTrt", all.x=TRUE)

#use ifelse statement to fill in the na values in the original carb fields
JarChem$JarpCO2<- ifelse((is.na(JarChem$JarpCO2)), paste(JarChem$EstJarpCO2), JarChem$JarpCO2)
JarChem$JarAlkCalc<- ifelse((is.na(JarChem$JarAlkCalc)), paste(JarChem$EstJarAlkCalc), JarChem$JarAlkCalc)
JarChem$JarHCO3<- ifelse((is.na(JarChem$JarHCO3)), paste(JarChem$EstJarHCO3), JarChem$JarHCO3)
JarChem$JarCO3<- ifelse((is.na(JarChem$JarCO3)), paste(JarChem$EstJarCO3), JarChem$JarCO3)
JarChem$JarCO2<- ifelse((is.na(JarChem$JarCO2)), paste(JarChem$EstJarCO2), JarChem$JarCO2)
JarChem$JarOmegaCalcite<- ifelse((is.na(JarChem$JarOmegaCalcite)), paste(JarChem$EstJarOmegaCalcite), JarChem$JarOmegaCalcite)
JarChem$JarOmegaArg<- ifelse((is.na(JarChem$JarOmegaArg)), paste(JarChem$EstJarOmegaArg), JarChem$JarOmegaArg)
#get the carb parameters to numeric
JarChem$JarpCO2<- as.numeric(JarChem$JarpCO2)
JarChem$JarAlkCalc<- as.numeric(JarChem$JarAlkCalc)
JarChem$JarHCO3<- as.numeric(JarChem$JarHCO3)
JarChem$JarCO3<- as.numeric(JarChem$JarCO3)
JarChem$JarCO2<- as.numeric(JarChem$JarCO2)
JarChem$JarOmegaCalcite<- as.numeric(JarChem$JarOmegaCalcite)
JarChem$JarOmegaArg<- as.numeric(JarChem$JarOmegaArg)

Look at the water chemistry data for adult tanks

#subset the tank data to only include B1 
Tank_WCB1<- subset(Tank_WC, Block == "B1")
Tank_WCB1$Tank<-droplevels(Tank_WCB1)$Tank #drop unused levels
#test to see if tanks are significantly different. 
boxplot(WCa_pHDIC~Tank, data=Tank_WCB1, na.omit=TRUE)

#check for significant differences in saturation state between tanks within treatments
Ca1<- lmer(WCa_pHDIC~ParentTrt + (1|Tank), data=Tank_WCB1)
## boundary (singular) fit: see ?isSingular
Ca2<- lmer(WCa_pHDIC~ 1 + (1|Tank), data= Tank_WCB1)
anova(Ca1, Ca2) #select Ca1
## refitting model(s) with ML (instead of REML)
## Data: Tank_WCB1
## Models:
## Ca2: WCa_pHDIC ~ 1 + (1 | Tank)
## Ca1: WCa_pHDIC ~ ParentTrt + (1 | Tank)
##     npar      AIC     BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## Ca2    3   9.9301  14.327 -1.9651   3.9301                         
## Ca1    4 -21.1940 -15.331 14.5970 -29.1940 33.124  1  8.646e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Ca3<- lm(WCa_pHDIC~ParentTrt, data=Tank_WCB1)
anova(Ca1, Ca3) #select Ca3
## refitting model(s) with ML (instead of REML)
## Data: Tank_WCB1
## Models:
## Ca3: WCa_pHDIC ~ ParentTrt
## Ca1: WCa_pHDIC ~ ParentTrt + (1 | Tank)
##     npar     AIC     BIC logLik deviance Chisq Df Pr(>Chisq)
## Ca3    3 -23.194 -18.797 14.597  -29.194                    
## Ca1    4 -21.194 -15.331 14.597  -29.194     0  1          1
#check assumptions
par(mfcol=c(2,2))
plot(Ca3) #seems to meet assumptions of homoscedasticity and linearity. Normality could look better. 

par(mfcol=c(1,1))
acf(resid(Ca3)) #this looks good. 

summary(Ca3) #There is a significant effect of ParentTrt, but not tank according to best model. 
## 
## Call:
## lm(formula = WCa_pHDIC ~ ParentTrt, data = Tank_WCB1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38809 -0.04330 -0.00918  0.06594  0.25937 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.20034    0.03959   30.32  < 2e-16 ***
## ParentTrt2800 -0.85889    0.05599  -15.34 9.66e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1584 on 30 degrees of freedom
## Multiple R-squared:  0.8869, Adjusted R-squared:  0.8832 
## F-statistic: 235.3 on 1 and 30 DF,  p-value: 9.659e-16
#try a transformation to see if it helps normality
Ca1log<- lmer(log(WCa_pHDIC)~ParentTrt+(1|Tank), data=Tank_WCB1)
Ca3log<- lm(log(WCa_pHDIC)~ParentTrt, data=Tank_WCB1)
anova(Ca1log, Ca3log) #still choose Ca3log
## refitting model(s) with ML (instead of REML)
## Data: Tank_WCB1
## Models:
## Ca3log: log(WCa_pHDIC) ~ ParentTrt
## Ca1log: log(WCa_pHDIC) ~ ParentTrt + (1 | Tank)
##        npar     AIC      BIC logLik deviance Chisq Df Pr(>Chisq)
## Ca3log    3 -15.952 -11.5545 10.976  -21.952                    
## Ca1log    4 -13.952  -8.0888 10.976  -21.952     0  1          1
par(mfcol=c(2,2))
plot(Ca3log) #now meets assumptions 

par(mfcol=c(1,1))
acf(resid(Ca3log))#looks good

summary(Ca3log)
## 
## Call:
## lm(formula = log(WCa_pHDIC) ~ ParentTrt, data = Tank_WCB1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.37383 -0.08255  0.00641  0.14423  0.29896 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.16588    0.04434   3.741 0.000773 ***
## ParentTrt2800 -1.25239    0.06270 -19.974  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1773 on 30 degrees of freedom
## Multiple R-squared:  0.9301, Adjusted R-squared:  0.9277 
## F-statistic:   399 on 1 and 30 DF,  p-value: < 2.2e-16
#create dataframe that has the mean values for each tank
TankMeans<- ddply(Tank_WCB1, .(Tank, Block, ParentTrt), numcolwise(mean, na.rm=T))

Join the data into dataframe for Egg analysis

#First dataframe will be for examining egg morphology for all eggs traced; multiple values for each adult
ParentInfo<- merge(x=Adult_Sample, y=Parent_ID, by="AdultID", all.x=TRUE)
ParentInfo<- merge(x=ParentInfo, y=TankMeans, by.x="TrtTankID", by.y="Tank", all.x=TRUE)
Eggs<- merge(x=Egg_Morphology, y=ParentInfo, by.x="FemaleID", by.y="ParentBlockID", all.x=TRUE)
Eggs$FemaleID<- as.factor(Eggs$FemaleID)
#subset the data to only include B1
EggsB1<- subset(Eggs, Block =="B1")

Visualize the data for eggs to start

#get the original range of the eggs
range(EggsB1$EggDiamum)
## [1] 15.35819 92.72958
#look at the egg measurements to see outliers
boxplot(EggDiamum~FemaleID, data=EggsB1)

boxplot(EggDiamum~ParentTrt, data=EggsB1)

Remove egg outliers

#check for outliers using code from https://www.r-bloggers.com/identify-describe-plot-and-remove-the-outliers-from-the-dataset/. The function uses Tukey's method to ID ouliers ranged above and below the 1.5*IQR. 
outlierKD <- function(dt, var) {
     var_name <- eval(substitute(var),eval(dt))
     na1 <- sum(is.na(var_name))
     m1 <- mean(var_name, na.rm = T)
     par(mfrow=c(2, 2), oma=c(0,0,3,0))
     boxplot(var_name, main="With outliers")
     hist(var_name, main="With outliers", xlab=NA, ylab=NA)
     outlier <- boxplot.stats(var_name)$out
     mo <- mean(outlier)
     var_name <- ifelse(var_name %in% outlier, NA, var_name)
     boxplot(var_name, main="Without outliers")
     hist(var_name, main="Without outliers", xlab=NA, ylab=NA)
     title("Outlier Check", outer=TRUE)
     na2 <- sum(is.na(var_name))
     cat("Outliers identified:", na2 - na1, "n")
     cat("Propotion (%) of outliers:", round((na2 - na1) / sum(!is.na(var_name))*100, 1), "n")
     cat("Mean of the outliers:", round(mo, 2), "n")
     m2 <- mean(var_name, na.rm = T)
     cat("Mean without removing outliers:", round(m1, 2), "n")
     cat("Mean if we remove outliers:", round(m2, 2), "n")
     response <- readline(prompt="Do you want to remove outliers and to replace with NA? [yes/no]: ")
     if(response == "y" | response == "yes"){
          dt[as.character(substitute(var))] <- invisible(var_name)
          assign(as.character(as.list(match.call())$dt), dt, envir = .GlobalEnv)
          cat("Outliers successfully removed", "n")
          return(invisible(dt))
     } else{
          dt[as.character(substitute(var))] <- invisible(var_name)
          assign(as.character(as.list(match.call())$dt), dt, envir = .GlobalEnv)
          cat("Outliers successfully removed", "n")
          return(invisible(dt))
     }
}
outlierKD(EggsB1, EggDiamum) #I opted to remove outliers for now. 

## Outliers identified: 48 nPropotion (%) of outliers: 11.5 nMean of the outliers: 57.63 nMean without removing outliers: 62.4 nMean if we remove outliers: 62.95 nDo you want to remove outliers and to replace with NA? [yes/no]: 
## Outliers successfully removed n
#Outliers identified: 48 nPropotion (%) of outliers: 11.5 nMean of the outliers: 57.63 nMean without removing outliers: 62.4 nMean if we remove outliers: 62.95 n
yes #to remove outliers from EggsB1
## Error in eval(expr, envir, enclos): object 'yes' not found
#count how many eggs are left per female after removing the outliers. 
outrem<- aggregate(EggDiamum~FemaleID, data=EggsB1, FUN=length)
#EF07_B1 had 10 eggs that were included the rest were removed. Keep EF07_B1 in the analysis

Visualize eggs without outliers

#Plot again without outliers
par(mfcol=c(1,1))
boxplot(EggDiamum~FemaleID, data=EggsB1) 

boxplot(EggDiamum~ParentTrt, data=EggsB1)

#get the range now
range(EggsB1$EggDiamum, na.rm=TRUE)
## [1] 49.95329 75.96955
#note that the eggs were filtered through a 70 um filter, but it appears that larger ones made it through.We are going to keep these in for now because they aren't that outrageously large. 
#plot histograms of the egg data with outliers removed
ggplot(EggsB1, aes(x=EggDiamum, color=ParentTrt))+
  geom_histogram(alpha=0.5, position="identity")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 48 rows containing non-finite values (stat_bin).

eggmeans<- aggregate(EggDiamum~ParentTrt, data=EggsB1, FUN=mean) #400 = 63.46233; 2800= 62.15902
eggmedians<- aggregate(EggDiamum~ParentTrt, data=EggsB1, FUN=median)#400 = 62.91200; 2800= 61.80109
view(eggmeans)
view(eggmedians)
#I will use the mean egg size. 

Make datasets without the egg size outliers

meanegg<- ddply(EggsB1, .(FemaleID), numcolwise(mean, na.rm=T)) #get mean egg size
#select only columns that we need for egg morph
meanegg<- subset(meanegg, select=c("FemaleID", "EggDiamum"))
FemAdults<- merge(x=ParentInfo, y= meanegg, by.x="ParentBlockID", by.y= "FemaleID", all.y=TRUE, all.x=TRUE)
#Next dataframe will be for one line per traced larva, includes all the data for that jar, including average egg size for that jar 
Larvae<-merge(x=Larvae_Morphology, y=JarChem, by="JarID", all.x=TRUE)
#use gather function in tidyr to make the Fert_QG Jar ID columns in long form vs. wide form
Fert_QG_long <- Fert_QG %>% tidyr::gather(Key, JarID, JarID1:JarID6)
#use the gathered data frame to merge with Larvae
Larvae<-merge(x=Larvae, y=Fert_QG_long, by.x= c("JarID", "CrossID"), by.y=c("JarID", "CrossID"), all.x=TRUE) #did the join by two columns so CrossID column wasn't duplicated
Larvae<- merge(x=Larvae, y=Cross, by="CrossID", all.x=TRUE)
Larvae<- merge(x=Larvae, y=FemAdults, by.x= "FemaleID", by.y="ParentBlockID", all.x=TRUE)
Larvae<- merge(x=Larvae, y= Larvae_Counts, by="JarID", all.x=TRUE)
#Remove jars without larvae
Larvae<- subset(Larvae, Larvae=="TRUE")
#get the growth per day for the larva
Larvae$GrowthPerDay<- (Larvae$LarvaeDiamum-Larvae$EggDiamum)/3
#Final dataframe will be one line per Jar, includes all data for the jar, including the average larva size per jar and cilia
LarvByJar<-merge(x=Larvae_Calcification, y=Larvae_Cilia, by="JarID", all.x=TRUE, all.y=TRUE)
LarvByJar<- merge(x=LarvByJar, y=JarChem, by="JarID", all.x=TRUE)
LarvByJar<- merge(x=LarvByJar, y=Fert_QG_long, by.x=c("JarID", "CrossID"), by.y=c("JarID","CrossID"), all.x=TRUE)
LarvByJar<- merge(x=LarvByJar, y=Cross, by="CrossID", all.x=TRUE)
LarvByJar<- merge(x=LarvByJar, y=FemAdults, by.x="FemaleID", by.y="ParentBlockID", all.x=TRUE)
LarvByJar<- merge(x=LarvByJar, y=Larvae_Counts, by="JarID", all.x=TRUE)
#get the means for larvae morphology and add to LarvByJar
MeanLarvMorph<- ddply(Larvae_Morphology, .(JarID), numcolwise(mean, na.rm=T))
SELarvMorph<- ddply(Larvae_Morphology, .(JarID), numcolwise(se, na.rm=T)) #get standard error of larva size
MeanLarvMorph<- full_join(MeanLarvMorph,SELarvMorph, by=c("JarID", "JarID"),suffix=c("","SE")) #join the mean and se larvae data
LarvByJar<- merge(x=LarvByJar, y=MeanLarvMorph, by="JarID", all.x=TRUE)
#Now let's add the survival data. See notes on "Larvae Survival" for more info on how we decided to count total larvae in the jar
#Make columns for the two larvae counts for 1_5 and 2
LarvByJar$v1_5SurvCount<- LarvByJar$SWRTotal+LarvByJar$F2LarvaeCount
LarvByJar$v2SurvCount<- LarvByJar$F1LarvaeCount+LarvByJar$F2LarvaeCount
#use ifelse statement to define the LarvaeSurvived column
LarvByJar$LarvaeSurvived<- ifelse(LarvByJar$ProtocolVersion =="1", paste(LarvByJar$TotalLarvae), ifelse(LarvByJar$ProtocolVersion =="1_5", paste(LarvByJar$v1_5SurvCount), paste(LarvByJar$v2SurvCount)))
LarvByJar$LarvaeSurvived<- as.numeric(LarvByJar$LarvaeSurvived)
## Warning: NAs introduced by coercion
#get the total larvae for the control jars and exposed jars for each family and then get the ratio of exposed to control to get an idea of how good a male and female pair are at larvae surviving. Can also use this to create reaction norms. 
SurvDat<- ddply(LarvByJar, .(CrossID, ParentTrt, JarTrt, MaleID, FemaleID, TrtTankID, AccTankID, Block), numcolwise(mean, na.rm=TRUE))
SESurvDat<- ddply(LarvByJar, .(CrossID, ParentTrt, JarTrt, Block), numcolwise(se, na.rm=TRUE)) #get the standard errors for each mean count
SESurvDat<- subset(SESurvDat, select= c(CrossID, JarTrt, LarvaeSurvived)) #rename the columns
colnames(SESurvDat)<- paste("SE", colnames(SESurvDat), sep="")
#merge the standard error data to the SurvDat dataframe
SurvDat<- merge(SurvDat, SESurvDat, by.x= c("CrossID", "JarTrt"), by.y= c("SECrossID", "SEJarTrt"))
#only keep the relevant fields because otherwise this is a huge dataframe
SurvDat<- subset(SurvDat, select= c(CrossID, JarTrt, ParentTrt, MaleID, FemaleID, TrtTankID, AccTankID, LarvaeSurvived, Block, SELarvaeSurvived, WCa_pHDIC)) 
SurvDat<- subset(SurvDat, Block=="B1")
SurvDat<- subset(SurvDat, CrossID!="")
#use the spread function to go from long form to wide form for the survival data
SurvDatToWide<- unite(SurvDat, Value, c(LarvaeSurvived, SELarvaeSurvived), sep="_")
SurvWideDat<- spread(SurvDatToWide, JarTrt, value= Value)
#now split the mean and SE apart again
SurvWideDat<- separate(SurvWideDat, Control, into= c("MeanLarvSurvivedCon", "SELarvSurvivedCon"), sep="_")
SurvWideDat<- separate(SurvWideDat, Exposed, into= c("MeanLarvSurvivedExp", "SELarvSurvivedExp"), sep="_")
#make the means and ses numeric
SurvWideDat$MeanLarvSurvivedExp<- as.numeric(SurvWideDat$MeanLarvSurvivedExp)
SurvWideDat$MeanLarvSurvivedCon<- as.numeric(SurvWideDat$MeanLarvSurvivedCon)
SurvWideDat$SELarvSurvivedExp<- as.numeric(SurvWideDat$SELarvSurvivedExp)
## Warning: NAs introduced by coercion
SurvWideDat$SELarvSurvivedCon<- as.numeric(SurvWideDat$SELarvSurvivedCon)
## Warning: NAs introduced by coercion
SurvWideDat$RatSurv<- SurvWideDat$MeanLarvSurvivedExp/SurvWideDat$MeanLarvSurvivedCon
SurvWideDat$SurvChange<- (SurvWideDat$MeanLarvSurvivedCon-SurvWideDat$MeanLarvSurvivedExp)/SurvWideDat$MeanLarvSurvivedCon
SurvWideDat$CrossID<- as.factor(SurvWideDat$CrossID)
SurvWideDat$FemaleID<- as.factor(SurvWideDat$FemaleID)
SurvWideDat$MaleID<- as.factor(SurvWideDat$MaleID)
LarvByJar<- subset(LarvByJar, Larvae=="TRUE")
#now remove the data we do not want to analyze.
LarvaeDat<- subset(Larvae, Block == "B1")
LarvaeDat$JarID<- as.factor(LarvaeDat$JarID)
LarvaeDat$AccTankID<- as.factor(LarvaeDat$AccTankID)
LarvaeDat$TrtTankID<- as.factor(LarvaeDat$TrtTankID)
LarvaeDat$JarSeatable<- as.factor(LarvaeDat$JarSeatable)
LarvaeDat$JarTrt<- as.factor(LarvaeDat$JarTrt)
LarvaeDat$ParentTrt<- as.factor(LarvaeDat$ParentTrt)
LarvaeDat$FemaleID<- as.factor(LarvaeDat$FemaleID)
LarvaeDat$MaleID<- as.factor(LarvaeDat$MaleID)
LarvaeDat$CrossID<- as.factor(LarvaeDat$CrossID)
LarvByJarDat<- subset(LarvByJar, Block == "B1")
LarvByJarDat$JarID<- as.factor(LarvByJarDat$JarID)
LarvByJarDat$AccTankID<- as.factor(LarvByJarDat$AccTankID)
LarvByJarDat$TrtTankID<- as.factor(LarvByJarDat$TrtTankID)
LarvByJarDat$JarSeatable<- as.factor(LarvByJarDat$JarSeatable)
LarvByJarDat$ParentTrt<- as.factor(LarvByJarDat$ParentTrt)
LarvByJarDat$FemaleID<- as.factor(LarvByJarDat$FemaleID)
LarvByJarDat$MaleID<- as.factor(LarvByJarDat$MaleID)
LarvByJarDat$CrossID<- as.factor(LarvByJarDat$CrossID)
LarvByJarDat$JarTrt<- as.factor(LarvByJarDat$JarTrt)
#make a column for the ratio of major to minor axis
LarvaeDat$MajMinRat<- LarvaeDat$LarvaeMajorAxisLengthPix/LarvaeDat$LarvaeMinorAxisLengthPix
LarvaeDat$PerimDiamRat<-LarvaeDat$LarvaePerimeterum/LarvaeDat$LarvaeDiamum

Analyze jar chemistry data

CalJar<- subset(LarvByJarDat, JarOmegaCalcite >0)
CalJar<- subset(CalJar, BottleChemistry=="TRUE" )#only use those that had their chem measured
#look at the saturation state for jars
boxplot(JarOmegaCalcite~JarTrt*JarSeatable, data=CalJar)

#I was having trouble with making this model: JarOmegaCalcite~JarTrt + (1|ParentTrt)+ (1|TimeFilter), so I looked at TimeFilter because I think that is the problem. 
Jar1<- lm(JarOmegaCalcite~TimeFilter, data=CalJar)
par(mfcol=c(2,2))
plot(Jar1)
## Warning: not plotting observations with leverage one:
##   1, 2, 3, 4, 5, 6, 7, 8
## Error in text.default(x, y, labels.id[ind], cex = cex.id, xpd = TRUE, : invalid 'pos' value
par(mfcol=c(1,1))

#looks like we can't use timefilter because it is too specific to each jar. Get the error message "observations with leverage one"; try using seatable instead
Jar2<- lmer(JarOmegaCalcite~JarTrt + (1|JarSeatable), data=CalJar)
Jar2step<- step(Jar2)
## Error in as_lmerModLmerTest(model): model not of class 'lmerMod': cannot coerce to class 'lmerModLmerTest
#get error model not of class lmerMod. I looked up this issue and turns out that error results when all random effects are dropped. See https://stackoverflow.com/questions/55638476/how-to-do-stepwise-model-with-random-effect-lme4-lmertest for more info if you want. 
JarFin<- lm(JarOmegaCalcite~JarTrt, data=CalJar)
#check assumptions 
par(mfcol=c(2,2))
plot(JarFin)
par(mfcol=c(1,1))
acf(resid(JarFin))
summary(JarFin)
## 
## Call:
## lm(formula = JarOmegaCalcite ~ JarTrt, data = CalJar)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.11383 -0.02406 -0.01337  0.02975  0.19438 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.0863     0.0413   26.30 4.69e-09 ***
## JarTrtExposed  -0.8486     0.0584  -14.53 4.93e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09234 on 8 degrees of freedom
## Multiple R-squared:  0.9635, Adjusted R-squared:  0.9589 
## F-statistic: 211.1 on 1 and 8 DF,  p-value: 4.933e-07
#Conclusion: JarSeatable doesn't need to be accounted for as a random effect in terms of water chem. 

Test for significant effect of parent treatment on egg size. Start with measured water chem

CalEggs<- subset(EggsB1, WCa_pHDIC != "NA")
CalEggs$TrtTankID<-droplevels(CalEggs)$TrtTankID
#start with models that have the response variable EggDiamum
eggdiam1<- lmer(EggDiamum~WCa_pHDIC*AdultLength + (1|FemaleID) + (1|TrtTankID) + (1|AccTankID), data=CalEggs)
## boundary (singular) fit: see ?isSingular
eggstep<-step(eggdiam1)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(eggstep)
## Backward reduced random-effect table:
## 
##                 Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                        8 -1210.6 2437.2                         
## (1 | TrtTankID)          1    7 -1210.6 2435.2  0.000  1          1    
## (1 | AccTankID)          2    6 -1210.6 2433.2  0.000  1          1    
## (1 | FemaleID)           0    5 -1236.3 2482.7 51.425  1  7.439e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                       Eliminated  Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## WCa_pHDIC:AdultLength          1  3.5785  3.5785     1 11.920  0.1926 0.6686
## AdultLength                    2  9.2137  9.2137     1 12.932  0.4958 0.4938
## WCa_pHDIC                      3 23.0328 23.0328     1 14.281  1.2394 0.2840
## 
## Model found:
## EggDiamum ~ (1 | FemaleID)
#final model chosen: EggDiamum~ 1 + (1|FemaleID)
eggdiamfin<- lmer(EggDiamum~1+(1|FemaleID), data=CalEggs)
#Check assumptions
par(mfcol=c(1,1))
qqnorm(resid(eggdiamfin))
qqline(resid(eggdiamfin))

plot(eggdiamfin)

acf(resid(eggdiamfin))

#the above looks fine, seems to meet all assumptions. 
summary(eggdiamfin)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EggDiamum ~ 1 + (1 | FemaleID)
##    Data: CalEggs
## 
## REML criterion at convergence: 2431.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8207 -0.6329 -0.1567  0.5843  3.4110 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  FemaleID (Intercept)  4.331   2.081   
##  Residual             18.586   4.311   
## Number of obs: 417, groups:  FemaleID, 16
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  62.8803     0.5638 15.1316   111.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Conclusion: Female matters but not treatment or adult length. Acclimation tank and treatment tank also don't matter.  
plot(EggDiamum~FemaleID, data=CalEggs) #we need to account for FemaleID, but not egg size in Larvae models. Can also make parallel models one with EggDiamum as a covariate and one with Female ID as a random effect

eggdiam1f<- lmer(EggDiamum~ParentTrt*AdultLength + (1|FemaleID) + (1|TrtTankID)+ (1|AccTankID), data=CalEggs)
## boundary (singular) fit: see ?isSingular
eggstepf<-step(eggdiam1f)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(eggstepf)
## Backward reduced random-effect table:
## 
##                 Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                        8 -1211.0 2438.1                         
## (1 | TrtTankID)          1    7 -1211.0 2436.1  0.000  1          1    
## (1 | AccTankID)          2    6 -1211.0 2434.1  0.000  1          1    
## (1 | FemaleID)           0    5 -1236.9 2483.8 51.771  1  6.236e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                       Eliminated  Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## ParentTrt:AdultLength          1  4.6817  4.6817     1 11.912  0.2519 0.6249
## AdultLength                    2  9.0581  9.0581     1 12.927  0.4874 0.4974
## ParentTrt                      3 20.2635 20.2635     1 14.231  1.0904 0.3138
## 
## Model found:
## EggDiamum ~ (1 | FemaleID)
#final model chosen: EggDiamum~ 1 + (1|FemaleID), which is the same as the eggdiamfin model
summary(eggdiamfin)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EggDiamum ~ 1 + (1 | FemaleID)
##    Data: CalEggs
## 
## REML criterion at convergence: 2431.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8207 -0.6329 -0.1567  0.5843  3.4110 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  FemaleID (Intercept)  4.331   2.081   
##  Residual             18.586   4.311   
## Number of obs: 417, groups:  FemaleID, 16
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  62.8803     0.5638 15.1316   111.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Choosing mean or median for larvae length data

ggplot(LarvaeDat, aes(x=LarvaeDiamum, color=JarTrt))+
  geom_histogram(alpha=0.5, position="identity") +
  facet_grid(~ParentTrt)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

larvmeans<- aggregate(LarvaeDiamum~ParentTrt*JarTrt, data=LarvaeDat, FUN=mean)#CP/CL : 75.24678; EP/CL: 75.46671; CP/EL: 64.22775; EP/EL: 65.62240
larvmedians<- aggregate(LarvaeDiamum~ParentTrt*JarTrt, data=LarvaeDat, FUN=median)#CP/CL : 75.29579; EP/CL: 75.52886; CP/EL: 64.27797; EP/EL: 65.82983
view(larvmeans)
view(larvmedians)
#I think we can use the mean for these data 

Look at the larvae diameter data

boxplot(LarvaeDiamum~FemaleID, data=LarvaeDat) #I looked specifically for weird things with EF07_B1 because that individual only had 10 eggs after removing outliers. But all seems fine! 

#create most complex model and then use the step function to see which model is best fit for the data; these are for ParentTrt and JarTrt as covariates
larvlen<-lmer(LarvaeDiamum~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarID)+(1|JarSeatable)+(1|AccTankID), data=LarvaeDat)
steppedlarvlen<-step(larvlen)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00444468 (tol = 0.002, component 1)
print(steppedlarvlen)
## Backward reduced random-effect table:
## 
##                   Eliminated npar  logLik   AIC     LRT Df Pr(>Chisq)    
## <none>                         11 -6759.2 13540                          
## (1 | AccTankID)            1   10 -6759.5 13539   0.698  1  0.4034678    
## (1 | CrossID)              2    9 -6760.6 13539   2.065  1  0.1507047    
## (1 | FemaleID)             0    8 -6777.8 13572  34.482  1  4.301e-09 ***
## (1 | MaleID)               0    8 -6766.8 13550  12.376  1  0.0004349 ***
## (1 | JarID)                0    8 -6856.8 13730 192.380  1  < 2.2e-16 ***
## (1 | JarSeatable)          0    8 -6763.7 13543   6.191  1  0.0128429 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated Sum Sq Mean Sq NumDF  DenDF F value
## WCa_pHDIC:JarOmegaCalcite          0 86.778  86.778     1 240.67   11.23
##                              Pr(>F)    
## WCa_pHDIC:JarOmegaCalcite 0.0009345 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## LarvaeDiamum ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) + 
##     (1 | MaleID) + (1 | JarID) + (1 | JarSeatable) + WCa_pHDIC:JarOmegaCalcite
#final model chosen by the lme4 step function: y~WCa_pHDIC*JarOmegaCalcite +(1|FemaleID)+(1|MaleID)+(1|JarID)+ (1|JarSeatable)
larvlenfin<- lmer(LarvaeDiamum~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|JarID)+ (1|JarSeatable), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvlenfin) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvlenfin))
qqline(resid(larvlenfin)) #has relatively long tails, doesn't seem to meet assumption of normality. Check on this

acf(resid(larvlenfin)) #this looks good 

summary(larvlenfin)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeDiamum ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +  
##     (1 | MaleID) + (1 | JarID) + (1 | JarSeatable)
##    Data: LarvaeDat
## 
## REML criterion at convergence: 13521.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8987 -0.5441  0.0495  0.6065  4.8949 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  JarID       (Intercept) 1.6344   1.2784  
##  FemaleID    (Intercept) 0.7649   0.8746  
##  MaleID      (Intercept) 0.3448   0.5872  
##  JarSeatable (Intercept) 0.1350   0.3675  
##  Residual                7.7275   2.7798  
## Number of obs: 2698, groups:  
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
## 
## Fixed effects:
##                           Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                63.7296     0.7644  26.2756  83.369  < 2e-16 ***
## WCa_pHDIC                  -2.1627     0.8158  24.3372  -2.651 0.013893 *  
## JarOmegaCalcite            10.8411     0.4941 240.7074  21.940  < 2e-16 ***
## WCa_pHDIC:JarOmegaCalcite   1.8102     0.5402 240.6723   3.351 0.000935 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC   -0.854                
## JarOmegClct -0.428  0.391         
## WC_HDIC:JOC  0.381 -0.437   -0.893
#Conclusion: signficant main effect of Parent and Jar calcite sat, and sig. interaction between parent and Jar saturation states

Use fixed effect for Jar treatment and adult tank as covariate

larvlenJF<-lmer(LarvaeDiamum~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarID)+(1|JarSeatable)+(1|AccTankID), data=LarvaeDat)
steppedlarvlenJF<-step(larvlenJF)
print(steppedlarvlenJF)
## Backward reduced random-effect table:
## 
##                   Eliminated npar  logLik   AIC     LRT Df Pr(>Chisq)    
## <none>                         11 -6757.4 13537                          
## (1 | AccTankID)            1   10 -6757.7 13535   0.702  1   0.402062    
## (1 | CrossID)              2    9 -6758.5 13535   1.634  1   0.201178    
## (1 | FemaleID)             0    8 -6776.1 13568  35.051  1  3.212e-09 ***
## (1 | MaleID)               0    8 -6766.1 13548  15.223  1  9.555e-05 ***
## (1 | JarID)                0    8 -6850.2 13716 183.331  1  < 2.2e-16 ***
## (1 | JarSeatable)          0    8 -6762.7 13541   8.292  1   0.003982 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## WCa_pHDIC:JarTrt          0 84.857  84.857     1 240.75  10.981 0.001062 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## LarvaeDiamum ~ WCa_pHDIC + JarTrt + (1 | FemaleID) + (1 | MaleID) + 
##     (1 | JarID) + (1 | JarSeatable) + WCa_pHDIC:JarTrt
#final model chosen by the lme4 step function: y~WCa_pHDIC*JarTrt +(1|FemaleID)+(1|MaleID)+(1|JarID)+ (1|JarSeatable)
larvlenfinJF<- lmer(LarvaeDiamum~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarID)+ (1|JarSeatable), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvlenfinJF) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvlenfinJF))
qqline(resid(larvlenfinJF)) #has relatively long tails, I think this should be acceptable

acf(resid(larvlenfinJF)) #this looks good 

summary(larvlenfinJF)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeDiamum ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | MaleID) +  
##     (1 | JarID) + (1 | JarSeatable)
##    Data: LarvaeDat
## 
## REML criterion at convergence: 13517.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9046 -0.5422  0.0507  0.6106  4.9021 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  JarID       (Intercept) 1.5778   1.2561  
##  FemaleID    (Intercept) 0.7587   0.8711  
##  MaleID      (Intercept) 0.3920   0.6261  
##  JarSeatable (Intercept) 0.1665   0.4081  
##  Residual                7.7275   2.7798  
## Number of obs: 2698, groups:  
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)              75.5290     0.7392  20.9914 102.172  < 2e-16 ***
## WCa_pHDIC                -0.2155     0.7824  19.3222  -0.275  0.78593    
## JarTrtExposed            -9.2409     0.4149 240.7741 -22.273  < 2e-16 ***
## WCa_pHDIC:JarTrtExposed  -1.5027     0.4535 240.7464  -3.314  0.00106 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrTrtE
## WCa_pHDIC   -0.841                
## JarTrtExpsd -0.281  0.259         
## WC_HDIC:JTE  0.251 -0.292   -0.893

Use fixed effect for both adult and larvae treatments

larvlenFE<-lmer(LarvaeDiamum~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarID)+(1|JarSeatable)+(1|TrtTankID)+(1|AccTankID), data=LarvaeDat)
## boundary (singular) fit: see ?isSingular
steppedlarvlenFE<-step(larvlenFE)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00645175 (tol = 0.002, component 1)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(steppedlarvlenFE)
## Backward reduced random-effect table:
## 
##                   Eliminated npar  logLik   AIC     LRT Df Pr(>Chisq)    
## <none>                         12 -6758.3 13541                          
## (1 | TrtTankID)            1   11 -6758.3 13539   0.000  1  0.9995923    
## (1 | AccTankID)            2   10 -6758.7 13537   0.641  1  0.4232503    
## (1 | CrossID)              3    9 -6759.4 13537   1.580  1  0.2087381    
## (1 | FemaleID)             0    8 -6776.6 13569  34.286  1  4.758e-09 ***
## (1 | MaleID)               0    8 -6767.0 13550  15.075  1  0.0001033 ***
## (1 | JarID)                0    8 -6852.0 13720 185.190  1  < 2.2e-16 ***
## (1 | JarSeatable)          0    8 -6763.6 13543   8.235  1  0.0041095 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## ParentTrt:JarTrt          0 75.391  75.391     1 240.74  9.7562 0.002006 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## LarvaeDiamum ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 | MaleID) + 
##     (1 | JarID) + (1 | JarSeatable) + ParentTrt:JarTrt
#final model chosen by the lme4 step function: y~ParentTrt*JarTrt +(1|FemaleID)+(1|MaleID)+(1|JarID)+ (1|JarSeatable)
larvlenfinFE<- lmer(LarvaeDiamum~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarID)+ (1|JarSeatable), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvlenfinFE) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvlenfinFE))
qqline(resid(larvlenfinFE)) #has relatively long tails, I think this should be acceptable

acf(resid(larvlenfinFE)) #this looks good 

summary(larvlenfinFE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeDiamum ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 | MaleID) +  
##     (1 | JarID) + (1 | JarSeatable)
##    Data: LarvaeDat
## 
## REML criterion at convergence: 13518.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9019 -0.5407  0.0506  0.6101  4.9087 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  JarID       (Intercept) 1.5895   1.2607  
##  FemaleID    (Intercept) 0.7497   0.8659  
##  MaleID      (Intercept) 0.3909   0.6252  
##  JarSeatable (Intercept) 0.1663   0.4078  
##  Residual                7.7275   2.7798  
## Number of obs: 2698, groups:  
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
## 
## Fixed effects:
##                             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                  75.2408     0.5109  15.6742 147.263  < 2e-16 ***
## ParentTrt2800                 0.2381     0.6448  19.2484   0.369  0.71593    
## JarTrtExposed               -11.0160     0.2564 240.7031 -42.959  < 2e-16 ***
## ParentTrt2800:JarTrtExposed   1.1720     0.3752 240.7394   3.123  0.00201 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.624               
## JarTrtExpsd -0.254  0.202        
## PrT2800:JTE  0.174 -0.293  -0.683

Visualize the larvae length data

#plot the data using the means for each jar
meanlarvjar<- ddply(LarvaeDat, .(JarID), numcolwise(mean, na.rm=T))
selarvjar<- ddply(LarvaeDat, .(JarID), numcolwise(se, na.rm=T))
ggplot(meanlarvjar,
       aes(fill = WCa_pHDIC, y = LarvaeDiamum, x = JarOmegaCalcite)) +
    geom_errorbar(aes(ymin= meanlarvjar$LarvaeDiamum-selarvjar$LarvaeDiamum, ymax=meanlarvjar$LarvaeDiamum+selarvjar$LarvaeDiamum),width=0.05, position=position_dodge(0))+
  geom_point(aes(colour=WCa_pHDIC)) +
  ggtitle("Larvae diameter by larvae jar saturation state")+ 
  ylab("Larvae Diameter (um)")+
  xlab("Larvae Jar Calcite Saturation State")+
  theme_classic()
## Warning: Use of `meanlarvjar$LarvaeDiamum` is discouraged. Use `LarvaeDiamum`
## instead.

## Warning: Use of `meanlarvjar$LarvaeDiamum` is discouraged. Use `LarvaeDiamum`
## instead.
## Warning: position_dodge requires non-overlapping x intervals

#the reason there is spread is that for the jars that we actually measured carbonate dynamics I used those values, for all other jars, I used the average of those measurements
#Try to visualize it in a more pleasant way
#get means and SE of larvae diameters
meanlarv<- ddply(LarvaeDat, .(JarTrt, ParentTrt), numcolwise(mean, na.rm=T))
selarv<- ddply(LarvaeDat, .(JarTrt, ParentTrt), numcolwise(se, na.rm=T))
#plot the mean Larvae Diameter +/- SE
ggplot(data = meanlarv, aes(x = JarOmegaCalcite, y = LarvaeDiamum, group=ParentTrt, fill=ParentTrt)) +
  geom_errorbar(aes(ymin= LarvaeDiamum-selarv$LarvaeDiamum, ymax=LarvaeDiamum+selarv$LarvaeDiamum),width=0.05, position=position_dodge(0))+
   geom_errorbarh(aes(xmin= JarOmegaCalcite-selarv$JarOmegaCalcite, xmax=JarOmegaCalcite+selarv$JarOmegaCalcite),width=0.05, position=position_dodge(0))+
  geom_point(aes(colour=ParentTrt))+
  labs(x = "Larvae Jar Saturation State (Calcite)", y = "Larvae Diameter (um)") +
  scale_color_manual(values = c("skyblue2", "red2"), name="Parent Treatment") +
theme_classic()
## Warning: Ignoring unknown parameters: width

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

#alternatively, I could visualize with bargraphs: 
#plot with parent treatment on the x-axis
ggplot(data = meanlarv, aes(x = ParentTrt , y = LarvaeDiamum, fill=JarTrt)) +
  geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
  labs(x = "Parent Treatment", y = "Larvae Diameter (um)") +
  scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA Exposed"), values = c("skyblue2", "red2")) +
  geom_errorbar(aes(ymin= meanlarv$LarvaeDiamum-selarv$LarvaeDiamum, ymax=meanlarv$LarvaeDiamum+selarv$LarvaeDiamum),width=0.2, position=position_dodge(.9))+
  theme_classic()
## Warning: Use of `meanlarv$LarvaeDiamum` is discouraged. Use `LarvaeDiamum`
## instead.
## Warning: Use of `meanlarv$LarvaeDiamum` is discouraged. Use `LarvaeDiamum`
## instead.

#now plot the difference between parental treatments for larvae diameter
crossdiff <- data.frame(tapply(LarvaeDat$LarvaeDiamum, list(LarvaeDat$CrossID, LarvaeDat$JarTrt), mean)) %>%
  rownames_to_column("CrossID") %>%
  mutate(CrossTemp = CrossID) %>%
  mutate(ParentTrt = ifelse(startsWith(CrossID, "E"), "Exposed", "Control")) %>%
  separate(CrossTemp, c("FemaleID", "MaleID", "BlockID", "Fourth"), sep="_") %>%
  mutate(Diff = Exposed - Control)
crossdiff$CrossID<- as.factor(crossdiff$CrossID)
crossdiff$FemaleID<-as.factor(crossdiff$FemaleID)
crossdiff$MaleID<-as.factor(crossdiff$MaleID)
#look at cross differences by Male and Female ID
ggplot(data = crossdiff, aes(x = MaleID, y = Diff)) + 
  geom_point(aes(colour=FemaleID)) + 
  labs(x = "Male ID", y = "Change in mean shell length (um) per family\nfrom control to OA conditions") + 
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(colour = "black"), panel.border=element_rect(colour="black", fill=NA))+
   theme(axis.text.x = element_text(angle = 90, vjust = 0.5))

#get the means of the differences to visualize overall differences
meandiff<- aggregate(Diff~ParentTrt, data=crossdiff, FUN=mean)
sediff<-  aggregate(Diff~ParentTrt, data=crossdiff, FUN=se)
ggplot(data = meandiff, aes(x = ParentTrt, y = Diff)) + 
  geom_bar(stat="identity", aes(fill="gray45"), position="dodge") + 
  labs(x = "Parental Treatment", y = "Change in mean shell length (um) per family\nfrom control to OA conditions") + 
  scale_x_discrete(breaks = c("Control", "Exposed"), labels = c("Control", "OA Exposed")) + 
  scale_fill_manual(values = "gray45") + 
  geom_errorbar(aes(ymin= Diff-sediff$Diff, ymax=Diff+sediff$Diff),width=0.2, position=position_dodge(.9))+
  theme(legend.position = "none", panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

diffttest<- t.test(Diff~ParentTrt, data=crossdiff, var.equal=T)
diffttest #mean shell length difference is significantly higher for OA parent treatment
## 
##  Two Sample t-test
## 
## data:  Diff by ParentTrt
## t = -2.6006, df = 43, p-value = 0.01271
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.1378060 -0.2703394
## sample estimates:
## mean in group Control mean in group Exposed 
##            -11.049847             -9.845775

Analyze at growth per day. This will likely be the main figure

boxplot(GrowthPerDay~ParentTrt*JarTrt, data=LarvaeDat)

larvgrow<-lmer(GrowthPerDay~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarID) +(1|JarSeatable)+(1|AccTankID) , data=LarvaeDat)
## boundary (singular) fit: see ?isSingular
steppedlarvgrow<- step(larvgrow)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
## Warning: Model failed to converge with 1 negative eigenvalue: -1.8e+00
print(steppedlarvgrow)
## Backward reduced random-effect table:
## 
##                   Eliminated npar  logLik    AIC     LRT Df Pr(>Chisq)    
## <none>                         11 -3809.8 7641.7                          
## (1 | AccTankID)            1   10 -3809.8 7639.7   0.000  1   0.999686    
## (1 | FemaleID)             0    9 -3828.6 7675.1  37.443  1  9.411e-10 ***
## (1 | MaleID)               0    9 -3811.7 7641.5   3.782  1   0.051815 .  
## (1 | CrossID)              0    9 -3814.3 7646.7   9.005  1   0.002692 ** 
## (1 | JarID)                0    9 -3891.5 7801.0 163.351  1  < 2.2e-16 ***
## (1 | JarSeatable)          0    9 -3813.1 7644.1   6.449  1   0.011100 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated Sum Sq Mean Sq NumDF  DenDF F value
## WCa_pHDIC:JarOmegaCalcite          0 10.242  10.242     1 221.49  11.928
##                              Pr(>F)    
## WCa_pHDIC:JarOmegaCalcite 0.0006624 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## GrowthPerDay ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) + 
##     (1 | MaleID) + (1 | CrossID) + (1 | JarID) + (1 | JarSeatable) + 
##     WCa_pHDIC:JarOmegaCalcite
#final model chosen by the lme4 step function: y~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|MaleID)+ (1|CrossID)+ (1|JarID) + (1|JarSeatable)
larvgrowfin<- lmer(GrowthPerDay~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|MaleID)+(1|CrossID) + (1|JarID) + (1|JarSeatable), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvgrowfin) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvgrowfin))
qqline(resid(larvgrowfin)) #has relatively long tails, based on scale I think it reasonably meets assumption of normality. 

acf(resid(larvgrowfin))

summary(larvgrowfin)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: GrowthPerDay ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +  
##     (1 | MaleID) + (1 | CrossID) + (1 | JarID) + (1 | JarSeatable)
##    Data: LarvaeDat
## 
## REML criterion at convergence: 7619.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9566 -0.5404  0.0513  0.6056  4.9473 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  JarID       (Intercept) 0.16985  0.4121  
##  CrossID     (Intercept) 0.02121  0.1456  
##  FemaleID    (Intercept) 0.43849  0.6622  
##  MaleID      (Intercept) 0.03583  0.1893  
##  JarSeatable (Intercept) 0.01485  0.1218  
##  Residual                0.85862  0.9266  
## Number of obs: 2698, groups:  
## JarID, 270; CrossID, 45; FemaleID, 15; MaleID, 11; JarSeatable, 3
## 
## Fixed effects:
##                           Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                 0.8459     0.4266  17.7762   1.983 0.063027 .  
## WCa_pHDIC                  -1.3679     0.4621  17.1680  -2.960 0.008705 ** 
## JarOmegaCalcite             3.6130     0.1611 221.3094  22.432  < 2e-16 ***
## WCa_pHDIC:JarOmegaCalcite   0.6082     0.1761 221.4933   3.454 0.000662 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC   -0.879                
## JarOmegClct -0.250  0.225         
## WC_HDIC:JOC  0.223 -0.251   -0.893
#Conclusion: Signficant main effect of tank and jar calcite saturation state and significant interaction.

Analyze larvae growth with adult trt as covariate and jar as fixed effect

larvgrowJF<-lmer(GrowthPerDay~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarID) +(1|JarSeatable)+(1|AccTankID) , data=LarvaeDat)
## boundary (singular) fit: see ?isSingular
steppedlarvgrowJF<- step(larvgrowJF)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
## Warning: Model failed to converge with 1 negative eigenvalue: -1.7e+00
print(steppedlarvgrowJF)
## Backward reduced random-effect table:
## 
##                   Eliminated npar  logLik    AIC     LRT Df Pr(>Chisq)    
## <none>                         11 -3808.0 7638.1                          
## (1 | AccTankID)            1   10 -3808.0 7636.1   0.000  1   1.000000    
## (1 | FemaleID)             0    9 -3827.0 7672.0  37.986  1  7.126e-10 ***
## (1 | MaleID)               0    9 -3810.6 7639.1   5.058  1   0.024518 *  
## (1 | CrossID)              0    9 -3812.3 7642.5   8.449  1   0.003653 ** 
## (1 | JarID)                0    9 -3886.6 7791.2 157.122  1  < 2.2e-16 ***
## (1 | JarSeatable)          0    9 -3812.3 7642.6   8.566  1   0.003426 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## WCa_pHDIC:JarTrt          0 9.9793  9.9793     1 221.43  11.623 0.0007741 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## GrowthPerDay ~ WCa_pHDIC + JarTrt + (1 | FemaleID) + (1 | MaleID) + 
##     (1 | CrossID) + (1 | JarID) + (1 | JarSeatable) + WCa_pHDIC:JarTrt
#final model chosen by the lme4 step function: y~WCa_pHDIC*JarTrt+(1|FemaleID)+(1|MaleID)+ (1|CrossID)+ (1|JarID) + (1|JarSeatable)
larvgrowfinJF<- lmer(GrowthPerDay~WCa_pHDIC*JarTrt+(1|FemaleID)+(1|MaleID)+(1|CrossID) + (1|JarID) + (1|JarSeatable), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvgrowfinJF) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvgrowfinJF))
qqline(resid(larvgrowfinJF)) #has relatively long tails, based on scale I think it reasonably meets assumption of normality. 

acf(resid(larvgrowfinJF))

summary(larvgrowfinJF)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: GrowthPerDay ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | MaleID) +  
##     (1 | CrossID) + (1 | JarID) + (1 | JarSeatable)
##    Data: LarvaeDat
## 
## REML criterion at convergence: 7616.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9576 -0.5394  0.0530  0.6100  4.9506 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  JarID       (Intercept) 0.16522  0.4065  
##  CrossID     (Intercept) 0.01806  0.1344  
##  FemaleID    (Intercept) 0.43793  0.6618  
##  MaleID      (Intercept) 0.04226  0.2056  
##  JarSeatable (Intercept) 0.01830  0.1353  
##  Residual                0.85862  0.9266  
## Number of obs: 2698, groups:  
## JarID, 270; CrossID, 45; FemaleID, 15; MaleID, 11; JarSeatable, 3
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)               4.7790     0.4222  17.0367  11.319 2.39e-09 ***
## WCa_pHDIC                -0.7149     0.4565  16.3546  -1.566 0.136483    
## JarTrtExposed            -3.0787     0.1356 221.2118 -22.703  < 2e-16 ***
## WCa_pHDIC:JarTrtExposed  -0.5054     0.1482 221.4272  -3.409 0.000774 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrTrtE
## WCa_pHDIC   -0.876                
## JarTrtExpsd -0.161  0.145         
## WC_HDIC:JTE  0.144 -0.163   -0.893

Analyze larvae growth with adult trt as covariate and jar as fixed effect

larvgrowFE<-lmer(GrowthPerDay~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarID) +(1|JarSeatable)+(1|TrtTankID)+(1|AccTankID) , data=LarvaeDat)
## boundary (singular) fit: see ?isSingular
steppedlarvgrowFE<- step(larvgrowFE)
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -9.5e+01
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -9.9e+00
## boundary (singular) fit: see ?isSingular
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.113534 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00209849 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.113534 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00235773 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.006054 (tol = 0.002, component 1)
print(steppedlarvgrowFE)
## Backward reduced random-effect table:
## 
##                   Eliminated npar  logLik    AIC     LRT Df Pr(>Chisq)    
## <none>                         12 -3809.2 7642.4                          
## (1 | AccTankID)            1   11 -3809.2 7640.4   0.000  1   1.000000    
## (1 | TrtTankID)            2   10 -3809.2 7638.4   0.000  1   0.999366    
## (1 | CrossID)              3    9 -3809.9 7637.7   1.310  1   0.252349    
## (1 | FemaleID)             0    8 -3905.5 7826.9 191.164  1  < 2.2e-16 ***
## (1 | MaleID)               0    8 -3818.2 7652.4  16.666  1  4.458e-05 ***
## (1 | JarID)                0    8 -3902.2 7820.4 184.710  1  < 2.2e-16 ***
## (1 | JarSeatable)          0    8 -3813.9 7643.9   8.156  1   0.004293 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## ParentTrt:JarTrt          0 8.3917  8.3917     1 241.78  9.7721 0.001989 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## GrowthPerDay ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 | MaleID) + 
##     (1 | JarID) + (1 | JarSeatable) + ParentTrt:JarTrt
#final model chosen by the lme4 step function: y~WCa_pHDIC*JarTrt+(1|FemaleID)+(1|MaleID)+ (1|CrossID)+ (1|JarID) + (1|JarSeatable)
larvgrowfinFE<- lmer(GrowthPerDay~ParentTrt*JarTrt+(1|FemaleID)+(1|MaleID)+(1|CrossID) + (1|JarID) + (1|JarSeatable), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvgrowfinFE) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvgrowfinFE))
qqline(resid(larvgrowfinFE)) #has relatively long tails, based on scale I think it reasonably meets assumption of normality. 

acf(resid(larvgrowfinFE))

summary(larvgrowfinFE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: GrowthPerDay ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 | MaleID) +  
##     (1 | CrossID) + (1 | JarID) + (1 | JarSeatable)
##    Data: LarvaeDat
## 
## REML criterion at convergence: 7618.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9541 -0.5394  0.0537  0.6073  4.9570 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  JarID       (Intercept) 0.16666  0.4082  
##  CrossID     (Intercept) 0.01776  0.1333  
##  FemaleID    (Intercept) 0.45139  0.6719  
##  MaleID      (Intercept) 0.04228  0.2056  
##  JarSeatable (Intercept) 0.01828  0.1352  
##  Residual                0.85861  0.9266  
## Number of obs: 2698, groups:  
## JarID, 270; CrossID, 45; FemaleID, 15; MaleID, 11; JarSeatable, 3
## 
## Fixed effects:
##                              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                   3.92983    0.27447  18.11960  14.318 2.54e-11 ***
## ParentTrt2800                 0.56798    0.38169  16.29530   1.488  0.15583    
## JarTrtExposed                -3.67574    0.08388 221.74326 -43.823  < 2e-16 ***
## ParentTrt2800:JarTrtExposed   0.39433    0.12271 221.44574   3.214  0.00151 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.661               
## JarTrtExpsd -0.155  0.111        
## PrT2800:JTE  0.106 -0.162  -0.684

Larvae growth per day visualization

#plot growth per day
ggplot(LarvaeDat,
       aes(fill = WCa_pHDIC, y = GrowthPerDay, x = JarOmegaCalcite)) +
  geom_point(aes(colour=WCa_pHDIC)) +
  ggtitle("Larvae growth by parent treatment")+ 
  ylab("Growth (um/day)")+
  theme_classic()

#get means and SE of larvae growth
meangrowth<- aggregate(GrowthPerDay~ParentTrt*JarTrt, data=LarvaeDat, FUN=mean)
segrowth<- aggregate(GrowthPerDay~ParentTrt*JarTrt, data=LarvaeDat, FUN=se)
meangrowth$JarTrt<- as.factor(meangrowth$JarTrt)
segrowth$JarTrt<- as.factor(segrowth$JarTrt)
#use a bargraph to plot the mean Larvae Diameter +/- SE
ggplot(data = meangrowth, aes(x = ParentTrt, y = GrowthPerDay, fill=JarTrt)) +
  geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
  labs(x = "Parent Treatment", y = "Larvae Growth (um/day)") +
  scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA Exposed"), values = c("skyblue2", "red2")) +
  geom_errorbar(aes(ymin= GrowthPerDay-segrowth$GrowthPerDay, ymax=GrowthPerDay+segrowth$GrowthPerDay),width=0.2, position=position_dodge(.9))+
  theme_classic()

#make plot like the Chirgwin et al. paper
par(mar = c(5, 4, 0.5, 0.5), bg = "transparent")
barCenters <- barplot(meangrowth$GrowthPerDay~meangrowth$JarTrt+meangrowth$ParentTrt,
                      space= c(0,.25), 
                      beside = TRUE, las = 1,
                      ylim = c(0, 5),
                      cex.names = 0.75,
                      ylab = expression(paste("Growth (", mu, "m/day)")),
                      xlab = "Larvae Treatment",
                      border = "black", axes = TRUE,
                      names=c("Control", "OA Exposed", "Control", "OA Exposed"),
                      legend.text = FALSE, col= c(rgb(0, 0, 1, alpha = 0.5), rgb(1, 0, 0, alpha = 0.5)))
arrows(x0=barCenters, y0=meangrowth$GrowthPerDay[c(1,3,2,4)] - segrowth$GrowthPerDay[c(1,3,2,4)], y1=meangrowth$GrowthPerDay[c(1,3,2,4)] + segrowth$GrowthPerDay[c(1,3,2,4)], angle=90, code=3)
rect(xleft = -0.2, 
     ybottom = 0, 
     xright = 2.37, 
     ytop = 5,
     border = TRUE,
     col = rgb(0, 0, 1, alpha = 0.15))
rect(xleft = 2.37, 
     ybottom = 0, 
     xright = 5.2, 
     ytop = 5,
     border = TRUE,
     col = rgb(1, 0, 0, alpha = 0.15))
barplot(meangrowth$GrowthPerDay~meangrowth$JarTrt+meangrowth$ParentTrt,
                      space= c(0,.25), 
                      beside = TRUE, las = 1,
                      ylim = c(0, 5),
                      cex.names = 0.75,
                      border = "black", axes = FALSE,
                      names=c("Control", "OA Exposed", "Control", "OA Exposed"),
                      legend.text = FALSE, col= c(rgb(0, 0, 1, alpha = 0.5), rgb(1, 0, 0, alpha = 0.5)), add=TRUE)

#now plot the difference between parental treatments for larvae growth
growthdiff <- data.frame(tapply(LarvaeDat$GrowthPerDay, list(LarvaeDat$CrossID, LarvaeDat$JarTrt), mean)) %>%
  rownames_to_column("CrossID") %>%
  mutate(CrossTemp = CrossID) %>%
  mutate(ParentTrt = ifelse(startsWith(CrossID, "E"), "Exposed", "Control")) %>%
  separate(CrossTemp, c("FemaleID", "MaleID", "BlockID", "Fourth"), sep="_") %>% mutate(Diff = Exposed - Control)
growthdiff$CrossID<- as.factor(growthdiff$CrossID)
growthdiff$FemaleID<-as.factor(growthdiff$FemaleID)
growthdiff$MaleID<-as.factor(growthdiff$MaleID)
#look at cross differences by Male and Female ID
ggplot(data = growthdiff, aes(x = MaleID, y = Diff)) + 
  geom_point(aes(colour=FemaleID)) + 
  labs(x = "Male ID", y = "Change in larvae growth (um/day) per family\nfrom control to OA conditions") + 
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(colour = "black"), panel.border=element_rect(colour="black", fill=NA))+
   theme(axis.text.x = element_text(angle = 90, vjust = 0.5))

#get the means of the differences to visualize overall differences
meangrdiff<- aggregate(Diff~ParentTrt, data=growthdiff, FUN=mean)
segrdiff<-  aggregate(Diff~ParentTrt, data=growthdiff, FUN=se)
ggplot(data = meangrdiff, aes(x = ParentTrt, y = Diff)) + 
  geom_bar(stat="identity", aes(fill="gray45"), position="dodge") + 
  labs(x = "Parental Treatment", y = "Change in larvae growth (um) per family\nfrom control to OA conditions") + 
  scale_x_discrete(breaks = c("Control", "Exposed"), labels = c("Control", "OA")) + 
  scale_fill_manual(values = "gray45") + 
  geom_errorbar(aes(ymin= Diff-segrdiff$Diff, ymax=Diff+segrdiff$Diff),width=0.2, position=position_dodge(.9))+
  theme(legend.position = "none", panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

diffgrttest<- t.test(Diff~ParentTrt, data=growthdiff, var.equal=T)
diffgrttest #significant effect of parent treatment on growth. 
## 
##  Two Sample t-test
## 
## data:  Diff by ParentTrt
## t = -2.6006, df = 43, p-value = 0.01271
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.71260201 -0.09011314
## sample estimates:
## mean in group Control mean in group Exposed 
##             -3.683282             -3.281925

Look at larvae area

boxplot(LarvaeAreaum2~ParentTrt*JarTrt, data=LarvaeDat)

#create most complex model and then use the step function to see which model is best fit for the data; these are for ParentTrt and JarTrt as covariates
larvarea<-lmer(LarvaeAreaum2~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) , data=LarvaeDat)
steppedlarvarea<-step(larvarea)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00443797 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00338512 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00332707 (tol = 0.002, component 1)
print(steppedlarvarea)
## Backward reduced random-effect table:
## 
##                   Eliminated npar logLik   AIC     LRT Df Pr(>Chisq)    
## <none>                         11 -18923 37868                          
## (1 | MaleID)               1   10 -18923 37867   1.097  1  0.2949993    
## (1 | AccTankID)            2    9 -18924 37866   1.464  1  0.2262316    
## (1 | JarSeatable)          3    8 -18925 37866   2.205  1  0.1375245    
## (1 | FemaleID)             0    7 -18928 37870   5.405  1  0.0200845 *  
## (1 | CrossID)              0    7 -18931 37876  12.013  1  0.0005284 ***
## (1 | JarID)                0    7 -19033 38079 214.833  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated Sum Sq Mean Sq NumDF  DenDF F value
## WCa_pHDIC:JarOmegaCalcite          0 376864  376864     1 223.37  5.9083
##                            Pr(>F)  
## WCa_pHDIC:JarOmegaCalcite 0.01586 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## LarvaeAreaum2 ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) + 
##     (1 | CrossID) + (1 | JarID) + WCa_pHDIC:JarOmegaCalcite
#final model chosen by the lme4 step function: y~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|CrossID)+(1|JarID)
larvareafin<- lmer(LarvaeAreaum2~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|CrossID)+(1|JarID), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvareafin) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvareafin))
qqline(resid(larvareafin)) #has relatively long tails, doesn't seem to meet assumption of normality. 

acf(resid(larvareafin)) #this looks fine to me. 

summary(larvareafin) #try transforming the data
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeAreaum2 ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +  
##     (1 | CrossID) + (1 | JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: 37850.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.5267 -0.5798  0.0408  0.6257  5.2994 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  JarID    (Intercept) 15333    123.82  
##  CrossID  (Intercept)  4909     70.06  
##  FemaleID (Intercept)  5197     72.09  
##  Residual             63786    252.56  
## Number of obs: 2698, groups:  JarID, 270; CrossID, 45; FemaleID, 15
## 
## Fixed effects:
##                           Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)                2660.72      60.11   24.02  44.267   <2e-16 ***
## WCa_pHDIC                  -151.83      65.62   23.94  -2.314   0.0296 *  
## JarOmegaCalcite            1142.11      46.93  223.22  24.337   <2e-16 ***
## WCa_pHDIC:JarOmegaCalcite   124.73      51.32  223.37   2.431   0.0159 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC   -0.893                
## JarOmegClct -0.517  0.462         
## WC_HDIC:JOC  0.461 -0.516   -0.893
larvareat<-lmer(sqrt(LarvaeAreaum2)~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) , data=LarvaeDat)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00224743 (tol = 0.002, component 1)
steppedlarvareat<-step(larvareat)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0020834 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00754904 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00754904 (tol = 0.002, component 1)
print(steppedlarvareat)
## Backward reduced random-effect table:
## 
##                   Eliminated npar  logLik   AIC     LRT Df Pr(>Chisq)    
## <none>                         11 -6182.9 12388                          
## (1 | AccTankID)            1   10 -6183.5 12387   1.118  1  0.2903515    
## (1 | MaleID)               2    9 -6184.4 12387   1.825  1  0.1767765    
## (1 | JarSeatable)          3    8 -6185.4 12387   2.114  1  0.1459162    
## (1 | FemaleID)             0    7 -6188.1 12390   5.308  1  0.0212333 *  
## (1 | CrossID)              0    7 -6191.3 12396  11.643  1  0.0006446 ***
## (1 | JarID)                0    7 -6286.0 12586 201.169  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated Sum Sq Mean Sq NumDF  DenDF F value
## WCa_pHDIC:JarOmegaCalcite          0 35.903  35.903     1 223.38  7.1837
##                             Pr(>F)   
## WCa_pHDIC:JarOmegaCalcite 0.007905 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## sqrt(LarvaeAreaum2) ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) + 
##     (1 | CrossID) + (1 | JarID) + WCa_pHDIC:JarOmegaCalcite
#final model chosen by the lme4 step function: y~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|CrossID)+(1|JarID)
larvareafint<- lmer(sqrt(LarvaeAreaum2)~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|CrossID)+(1|JarID), data=LarvaeDat)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00754904 (tol = 0.002, component 1)
#got warning message that model failed to converged. I read this stack overflow talking about checking to see if it is acceptable: https://stats.stackexchange.com/questions/97929/lmer-model-fails-to-converge
TestConverge <- with(larvareafint@optinfo$derivs,solve(Hessian,gradient))
max(abs(TestConverge))
## [1] 0.0004658432
#this value seems small enough that it's acceptable
plot(larvareafint) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvareafint))
qqline(resid(larvareafint)) #this seems to have improved it

acf(resid(larvareafint)) #this looks fine to me. 

summary(larvareafint)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sqrt(LarvaeAreaum2) ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +  
##     (1 | CrossID) + (1 | JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: 12370.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.5313 -0.5499  0.0578  0.6247  5.2230 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  JarID    (Intercept) 1.1442   1.0697  
##  CrossID  (Intercept) 0.3637   0.6031  
##  FemaleID (Intercept) 0.3846   0.6202  
##  Residual             4.9979   2.2356  
## Number of obs: 2698, groups:  JarID, 270; CrossID, 45; FemaleID, 15
## 
## Fixed effects:
##                           Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                51.7442     0.5194  24.2032  99.624  < 2e-16 ***
## WCa_pHDIC                  -1.4302     0.5670  24.1239  -2.522  0.01866 *  
## JarOmegaCalcite             9.8236     0.4084 223.2301  24.056  < 2e-16 ***
## WCa_pHDIC:JarOmegaCalcite   1.1968     0.4465 223.3805   2.680  0.00791 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC   -0.893                
## JarOmegClct -0.520  0.465         
## WC_HDIC:JOC  0.464 -0.519   -0.893
## convergence code: 0
## Model failed to converge with max|grad| = 0.00754904 (tol = 0.002, component 1)
#significant main effects of tank and jar treatments and significant interaction

Analyze surface area with Parent Trt as covariate and Jar Treatment as fixed effect

larvareatJF<-lmer(sqrt(LarvaeAreaum2)~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) , data=LarvaeDat)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00357282 (tol = 0.002, component 1)
steppedlarvareatJF<-step(larvareatJF)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00296653 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00461815 (tol = 0.002, component 1)
print(steppedlarvareatJF)
## Backward reduced random-effect table:
## 
##                   Eliminated npar  logLik   AIC     LRT Df Pr(>Chisq)    
## <none>                         11 -6181.1 12384                          
## (1 | AccTankID)            1   10 -6181.6 12383   1.111  1  0.2918610    
## (1 | CrossID)              2    9 -6182.9 12384   2.640  1  0.1042159    
## (1 | FemaleID)             0    8 -6201.7 12419  37.510  1  9.095e-10 ***
## (1 | MaleID)               0    8 -6190.0 12396  14.136  1  0.0001701 ***
## (1 | JarSeatable)          0    8 -6184.7 12385   3.444  1  0.0634889 .  
## (1 | JarID)                0    8 -6292.4 12601 218.856  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)  
## WCa_pHDIC:JarTrt          0  32.41   32.41     1 240.67  6.4848 0.0115 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## sqrt(LarvaeAreaum2) ~ WCa_pHDIC + JarTrt + (1 | FemaleID) + (1 | 
##     MaleID) + (1 | JarSeatable) + (1 | JarID) + WCa_pHDIC:JarTrt
#final model chosen by the lme4 step function: y~WCa_pHDIC*JarTrt+(1|FemaleID)+(1|MaleID)+(1|JarID)+(1|JarSeatable)
larvareafintJF<- lmer(sqrt(LarvaeAreaum2)~WCa_pHDIC*JarTrt+(1|FemaleID)+(1|MaleID)+(1|JarID)+(1|JarSeatable), data=LarvaeDat)
plot(larvareafintJF) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvareafintJF))
qqline(resid(larvareafintJF)) #this seems to have improved it

acf(resid(larvareafintJF)) #this looks fine to me. 

summary(larvareafintJF)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sqrt(LarvaeAreaum2) ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 |  
##     MaleID) + (1 | JarID) + (1 | JarSeatable)
##    Data: LarvaeDat
## 
## REML criterion at convergence: 12365.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4593 -0.5425  0.0603  0.6299  5.1752 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  JarID       (Intercept) 1.16250  1.0782  
##  FemaleID    (Intercept) 0.55911  0.7477  
##  MaleID      (Intercept) 0.26504  0.5148  
##  JarSeatable (Intercept) 0.06002  0.2450  
##  Residual                4.99788  2.2356  
## Number of obs: 2698, groups:  
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)              62.4349     0.6082  20.9581 102.663   <2e-16 ***
## WCa_pHDIC                -0.1523     0.6598  19.3614  -0.231   0.8199    
## JarTrtExposed            -8.3685     0.3489 240.6878 -23.986   <2e-16 ***
## WCa_pHDIC:JarTrtExposed  -0.9711     0.3813 240.6713  -2.547   0.0115 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrTrtE
## WCa_pHDIC   -0.863                
## JarTrtExpsd -0.287  0.259         
## WC_HDIC:JTE  0.257 -0.291   -0.893

Analyze surface area with Parent and Jar Treatment as fixed effects

larvareatFE<-lmer(sqrt(LarvaeAreaum2)~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) + (1|TrtTankID), data=LarvaeDat)
## boundary (singular) fit: see ?isSingular
steppedlarvareatFE<-step(larvareatFE)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00347465 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00448054 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00448054 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00238224 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0126724 (tol = 0.002, component 1)
print(steppedlarvareatFE)
## Backward reduced random-effect table:
## 
##                   Eliminated npar  logLik   AIC     LRT Df Pr(>Chisq)    
## <none>                         12 -6182.0 12388                          
## (1 | TrtTankID)            1   11 -6182.0 12386   0.000  1  0.9977038    
## (1 | AccTankID)            2   10 -6182.5 12385   1.085  1  0.2975686    
## (1 | CrossID)              3    9 -6183.8 12386   2.571  1  0.1088464    
## (1 | FemaleID)             0    8 -6202.1 12420  36.662  1  1.405e-09 ***
## (1 | MaleID)               0    8 -6190.8 12398  13.997  1  0.0001831 ***
## (1 | JarSeatable)          0    8 -6185.5 12387   3.418  1  0.0644882 .  
## (1 | JarID)                0    8 -6294.2 12604 220.780  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## ParentTrt:JarTrt          0 26.811  26.811     1 240.66  5.3644 0.02139 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## sqrt(LarvaeAreaum2) ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 | 
##     MaleID) + (1 | JarSeatable) + (1 | JarID) + ParentTrt:JarTrt
#final model chosen by the lme4 step function: y~ParentTrt*JarTrt+(1|FemaleID)+(1|MaleID)+(1|JarID)+(1|JarSeatable)
larvareafintFE<- lmer(sqrt(LarvaeAreaum2)~ParentTrt*JarTrt+(1|FemaleID)+(1|MaleID)+(1|JarID)+(1|JarSeatable), data=LarvaeDat)
plot(larvareafintFE) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvareafintFE))
qqline(resid(larvareafintFE)) #this seems to have improved it

acf(resid(larvareafintFE)) #this looks fine to me. 

summary(larvareafintFE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sqrt(LarvaeAreaum2) ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 |  
##     MaleID) + (1 | JarID) + (1 | JarSeatable)
##    Data: LarvaeDat
## 
## REML criterion at convergence: 12367.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4562 -0.5434  0.0596  0.6310  5.1810 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  JarID       (Intercept) 1.17016  1.0817  
##  FemaleID    (Intercept) 0.55120  0.7424  
##  MaleID      (Intercept) 0.26425  0.5141  
##  JarSeatable (Intercept) 0.05998  0.2449  
##  Residual                4.99788  2.2356  
## Number of obs: 2698, groups:  
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
## 
## Fixed effects:
##                             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                  62.2179     0.4069  18.4716 152.903   <2e-16 ***
## ParentTrt2800                 0.1959     0.5433  19.2917   0.361   0.7223    
## JarTrtExposed                -9.5031     0.2156 240.6370 -44.077   <2e-16 ***
## ParentTrt2800:JarTrtExposed   0.7307     0.3155 240.6577   2.316   0.0214 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.659               
## JarTrtExpsd -0.269  0.201        
## PrT2800:JTE  0.184 -0.292  -0.683

Larvae Area Visualization

#plot the data using the means for each jar
ggplot(LarvaeDat,
       aes(fill = WCa_pHDIC, y = LarvaeAreaum2, x = JarOmegaCalcite)) +
  geom_point(aes(shape= Block, colour=WCa_pHDIC)) +
  ggtitle("Larvae area by larvae treatment")+ 
  ylab("Larvae Area (um2)")+
  theme_classic()

#get means and SE of larvae areas
meanarea<- aggregate(LarvaeAreaum2~ParentTrt*JarTrt, data=LarvaeDat, FUN=mean)
searea<- aggregate(LarvaeAreaum2~ParentTrt*JarTrt, data=LarvaeDat, FUN=se)
#use a bargraph to plot the mean Larvae area +/- SE
ggplot(data = meanarea, aes(x = ParentTrt, y = LarvaeAreaum2, fill=JarTrt)) +
  geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
  labs(x = "Parent Treatment", y = "Larvae Area (um2)") +
  scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA"), values = c("skyblue2", "red2")) +
  geom_errorbar(aes(ymin= meanarea$LarvaeAreaum2-searea$LarvaeAreaum2, ymax=meanarea$LarvaeAreaum2+searea$LarvaeAreaum2),width=0.2, position=position_dodge(.9))+
  theme_classic()
## Warning: Use of `meanarea$LarvaeAreaum2` is discouraged. Use `LarvaeAreaum2`
## instead.

## Warning: Use of `meanarea$LarvaeAreaum2` is discouraged. Use `LarvaeAreaum2`
## instead.

#now plot the difference between parental treatments for larvae area
areadiff <- data.frame(tapply(LarvaeDat$LarvaeAreaum2, list(LarvaeDat$CrossID, LarvaeDat$JarTrt), mean)) %>%
  rownames_to_column("CrossID") %>%
  mutate(CrossTemp = CrossID) %>%
  mutate(ParentTrt = ifelse(startsWith(CrossID, "E"), "Exposed", "Control")) %>%
  separate(CrossTemp, c("FemaleID", "MaleID", "BlockID", "Fourth"), sep="_")%>%
  mutate(Diff = Exposed - Control)
areadiff$CrossID<- as.factor(areadiff$CrossID)
areadiff$FemaleID<-as.factor(areadiff$FemaleID)
areadiff$MaleID<-as.factor(areadiff$MaleID)
#look at cross differences by Male and Female ID
ggplot(data = areadiff, aes(x = MaleID, y = Diff)) + 
  geom_point(aes(colour=FemaleID)) + 
  labs(x = "Male ID", y = "Change in mean shell area (um2) per family\nfrom control to OA conditions") + 
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(colour = "black"), panel.border=element_rect(colour="black", fill=NA))+
   theme(axis.text.x = element_text(angle = 90, vjust = 0.5))

#get the means of the differences to visualize overall differences
meanareadiff<- aggregate(Diff~ParentTrt, data=areadiff, FUN=mean)
seareadiff<-  aggregate(Diff~ParentTrt, data=areadiff, FUN=se)
ggplot(data = meanareadiff, aes(x = ParentTrt, y = Diff)) + 
  geom_bar(stat="identity", aes(fill="gray45"), position="dodge") + 
  labs(x = "Parental Treatment", y = "Change in mean shell area (um2) per family\nfrom control to OA conditions") + 
  scale_x_discrete(breaks = c("Control", "Exposed"), labels = c("Control", "OA")) + 
  scale_fill_manual(values = "gray45") + 
  geom_errorbar(aes(ymin= meanareadiff$Diff-seareadiff$Diff, ymax=meanareadiff$Diff+seareadiff$Diff),width=0.2, position=position_dodge(.9))+
  theme(legend.position = "none", panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))
## Warning: Use of `meanareadiff$Diff` is discouraged. Use `Diff` instead.
## Warning: Use of `meanareadiff$Diff` is discouraged. Use `Diff` instead.

diffareattest<- t.test(Diff~ParentTrt, data=areadiff)
diffareattest #mean shell area difference is not significantly higher for OA parent treatment
## 
##  Welch Two Sample t-test
## 
## data:  Diff by ParentTrt
## t = -1.7272, df = 39.691, p-value = 0.09192
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -168.53796   13.23648
## sample estimates:
## mean in group Control mean in group Exposed 
##             -1093.243             -1015.592

Characterize the larvae by the ratio of major and minor axis

boxplot(MajMinRat~JarTrt*ParentTrt, data=LarvaeDat)

larvrat<-lmer(MajMinRat~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) +(1|CrossID), data=LarvaeDat)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00219325 (tol = 0.002, component 1)
larvratstep<- step(larvrat)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(larvratstep)
## Backward reduced random-effect table:
## 
##                   Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                         11 4750.1 -9478.2                         
## (1 | AccTankID)            1   10 4750.1 -9480.1  0.047  1     0.8285    
## (1 | JarSeatable)          2    9 4750.0 -9482.0  0.119  1     0.7298    
## (1 | MaleID)               3    8 4749.8 -9483.6  0.395  1     0.5299    
## (1 | CrossID)              4    7 4749.0 -9483.9  1.693  1     0.1932    
## (1 | FemaleID)             0    6 4738.3 -9464.6 21.363  1  3.801e-06 ***
## (1 | JarID)                0    6 4731.8 -9451.6 34.298  1  4.728e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated    Sum Sq   Mean Sq NumDF  DenDF F value
## WCa_pHDIC:JarOmegaCalcite          0 0.0068918 0.0068918     1 252.23  4.2751
##                            Pr(>F)  
## WCa_pHDIC:JarOmegaCalcite 0.03969 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## MajMinRat ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) + (1 | 
##     JarID) + WCa_pHDIC:JarOmegaCalcite
larvratfin<- lmer(MajMinRat~WCa_pHDIC*JarOmegaCalcite +(1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvratfin) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvratfin))
qqline(resid(larvratfin)) #has relatively long tails, doesn't seem to meet assumption of normality. try a transformation

acf(resid(larvratfin)) #this looks fine.

summary(larvratfin) #significant effect of parent trt and significant interaction between parent treatment and larval treatment
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MajMinRat ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) + (1 |  
##     JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: -9497.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6168 -0.5759 -0.0263  0.5463  4.7286 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  JarID    (Intercept) 1.078e-04 0.010382
##  FemaleID (Intercept) 5.034e-05 0.007095
##  Residual             1.612e-03 0.040151
## Number of obs: 2698, groups:  JarID, 270; FemaleID, 15
## 
## Fixed effects:
##                             Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                 1.171609   0.005785  31.004778 202.522   <2e-16 ***
## WCa_pHDIC                  -0.013115   0.006313  30.846657  -2.077   0.0462 *  
## JarOmegaCalcite            -0.008297   0.005224 252.351564  -1.588   0.1135    
## WCa_pHDIC:JarOmegaCalcite   0.011805   0.005710 252.231250   2.068   0.0397 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC   -0.893                
## JarOmegClct -0.598  0.534         
## WC_HDIC:JOC  0.533 -0.597   -0.893
#try transforming the data
larvratt<-lmer(MajMinRat^1/3~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) +(1|CrossID), data=LarvaeDat)
larvratstept<- step(larvratt)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00263927 (tol = 0.002, component 1)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(larvratstept)
## Backward reduced random-effect table:
## 
##                   Eliminated npar logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                         11 7709.8 -15398                         
## (1 | AccTankID)            1   10 7709.7 -15400  0.047  1     0.8285    
## (1 | JarSeatable)          2    9 7709.7 -15401  0.119  1     0.7298    
## (1 | MaleID)               3    8 7709.5 -15403  0.395  1     0.5299    
## (1 | CrossID)              4    7 7708.6 -15403  1.693  1     0.1932    
## (1 | FemaleID)             0    6 7697.9 -15384 21.363  1  3.801e-06 ***
## (1 | JarID)                0    6 7691.5 -15371 34.298  1  4.728e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated     Sum Sq    Mean Sq NumDF  DenDF F value
## WCa_pHDIC:JarOmegaCalcite          0 0.00076575 0.00076575     1 252.23  4.2751
##                            Pr(>F)  
## WCa_pHDIC:JarOmegaCalcite 0.03969 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## MajMinRat^1/3 ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) + 
##     (1 | JarID) + WCa_pHDIC:JarOmegaCalcite
larvratfint<- lmer((MajMinRat^1/3)~WCa_pHDIC*JarOmegaCalcite +(1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvratfint) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvratfint))
qqline(resid(larvratfint)) #has relatively long tails, doesn't seem to meet assumption of normality despite transformation

summary(larvratfint) #significant effect of parent trt and significant interaction between parent treatment and larval treatment
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: (MajMinRat^1/3) ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +  
##     (1 | JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: -15417.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6168 -0.5759 -0.0263  0.5463  4.7286 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  JarID    (Intercept) 1.198e-05 0.003461
##  FemaleID (Intercept) 5.593e-06 0.002365
##  Residual             1.791e-04 0.013384
## Number of obs: 2698, groups:  JarID, 270; FemaleID, 15
## 
## Fixed effects:
##                             Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                 0.390536   0.001928  31.004777 202.522   <2e-16 ***
## WCa_pHDIC                  -0.004372   0.002104  30.846656  -2.077   0.0462 *  
## JarOmegaCalcite            -0.002766   0.001741 252.351548  -1.588   0.1135    
## WCa_pHDIC:JarOmegaCalcite   0.003935   0.001903 252.231235   2.068   0.0397 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC   -0.893                
## JarOmegClct -0.598  0.534         
## WC_HDIC:JOC  0.533 -0.597   -0.893
#I have tried transforming with reciprocal, square root, square, and cube root none of them have helped with normality. No matter what though, still get a significant effect of parent treatment and significant interaction

Analyze the larvae by the ratio of major and minor axis with parent as a covariate and larvae as a fixed effect

larvratJF<-lmer(MajMinRat~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) +(1|CrossID), data=LarvaeDat)
larvratstepJF<- step(larvratJF)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(larvratstepJF)
## Backward reduced random-effect table:
## 
##                   Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                         11 4749.8 -9477.5                         
## (1 | AccTankID)            1   10 4749.7 -9479.5  0.046  1     0.8293    
## (1 | JarSeatable)          2    9 4749.7 -9481.4  0.114  1     0.7355    
## (1 | MaleID)               3    8 4749.5 -9483.0  0.397  1     0.5289    
## (1 | CrossID)              4    7 4748.6 -9483.3  1.705  1     0.1917    
## (1 | FemaleID)             0    6 4737.9 -9463.9 21.416  1  3.697e-06 ***
## (1 | JarID)                0    6 4731.5 -9451.0 34.285  1  4.760e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated    Sum Sq   Mean Sq NumDF  DenDF F value  Pr(>F)  
## WCa_pHDIC:JarTrt          0 0.0069002 0.0069002     1 252.22  4.2804 0.03957 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## MajMinRat ~ WCa_pHDIC + JarTrt + (1 | FemaleID) + (1 | JarID) + 
##     WCa_pHDIC:JarTrt
larvratfinJF<- lmer(MajMinRat~WCa_pHDIC*JarTrt +(1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvratfinJF) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvratfinJF))
qqline(resid(larvratfinJF)) #has relatively long tails, doesn't seem to meet assumption of normality. try a transformation

acf(resid(larvratfinJF)) #this looks fine.

summary(larvratfinJF)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MajMinRat ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: -9497.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6172 -0.5756 -0.0274  0.5463  4.7283 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  JarID    (Intercept) 1.078e-04 0.010381
##  FemaleID (Intercept) 5.043e-05 0.007101
##  Residual             1.612e-03 0.040151
## Number of obs: 2698, groups:  JarID, 270; FemaleID, 15
## 
## Fixed effects:
##                           Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)              1.163e+00  5.144e-03  1.958e+01 226.008   <2e-16 ***
## WCa_pHDIC               -2.818e-04  5.627e-03  1.966e+01  -0.050   0.9606    
## JarTrtExposed            7.059e-03  4.439e-03  2.523e+02   1.590   0.1130    
## WCa_pHDIC:JarTrtExposed -1.003e-02  4.850e-03  2.522e+02  -2.069   0.0396 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrTrtE
## WCa_pHDIC   -0.893                
## JarTrtExpsd -0.432  0.386         
## WC_HDIC:JTE  0.386 -0.434   -0.893

Analyze the larvae by the ratio of major and minor axis with parent and larvae as fixed effects

larvratFE<-lmer(MajMinRat~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|TrtTankID)+(1|AccTankID) +(1|CrossID), data=LarvaeDat)
larvratstepFE<- step(larvratFE)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -2.9e+01
## boundary (singular) fit: see ?isSingular
print(larvratstepFE)
## Backward reduced random-effect table:
## 
##                   Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                         12 4749.5 -9475.1                         
## (1 | TrtTankID)            1   11 4749.5 -9477.1  0.000  1     0.9991    
## (1 | AccTankID)            2   10 4749.5 -9479.0  0.047  1     0.8282    
## (1 | JarSeatable)          3    9 4749.5 -9480.9  0.116  1     0.7331    
## (1 | MaleID)               4    8 4749.3 -9482.5  0.392  1     0.5312    
## (1 | CrossID)              5    7 4748.4 -9482.8  1.730  1     0.1884    
## (1 | FemaleID)             0    6 4737.4 -9462.8 22.000  1  2.727e-06 ***
## (1 | JarID)                0    6 4731.4 -9450.8 34.009  1  5.486e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated    Sum Sq   Mean Sq NumDF  DenDF F value  Pr(>F)  
## ParentTrt:JarTrt          0 0.0076797 0.0076797     1 252.22  4.7639 0.02999 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## MajMinRat ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 | JarID) + 
##     ParentTrt:JarTrt
larvratfinFE<- lmer(MajMinRat~ParentTrt*JarTrt +(1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvratfinFE) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvratfinFE))
qqline(resid(larvratfinFE)) #has relatively long tails, doesn't seem to meet assumption of normality

acf(resid(larvratfinFE)) #this looks fine.

summary(larvratfinFE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MajMinRat ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 | JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: -9496.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6210 -0.5758 -0.0245  0.5460  4.7281 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  JarID    (Intercept) 0.0001072 0.010356
##  FemaleID (Intercept) 0.0000514 0.007169
##  Residual             0.0016121 0.040151
## Number of obs: 2698, groups:  JarID, 270; FemaleID, 15
## 
## Fixed effects:
##                               Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                  1.163e+00  3.195e-03  1.966e+01 363.863   <2e-16
## ParentTrt2800               -3.787e-04  4.670e-03  1.955e+01  -0.081   0.9362
## JarTrtExposed               -5.217e-03  2.733e-03  2.521e+02  -1.909   0.0574
## ParentTrt2800:JarTrtExposed  8.730e-03  4.000e-03  2.522e+02   2.183   0.0300
##                                
## (Intercept)                 ***
## ParentTrt2800                  
## JarTrtExposed               .  
## ParentTrt2800:JarTrtExposed *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.684               
## JarTrtExpsd -0.434  0.297        
## PrT2800:JTE  0.296 -0.431  -0.683

Visualize the ratio of the axes

#get the mean ratio for the treatments to plot
meanrat<- ddply(LarvaeDat, .(JarTrt, ParentTrt), numcolwise(mean, na.rm=T))
serat<- ddply(LarvaeDat, .(JarTrt, ParentTrt), numcolwise(se, na.rm=T))
ggplot(LarvaeDat, aes(x=MajMinRat, color=JarTrt))+
  geom_histogram(alpha=0.5, position="identity") +
  facet_grid(~ParentTrt)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(data = meanrat, aes(x = ParentTrt, y = MajMinRat, fill=JarTrt)) +
  geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
  labs(x = "Parent Treatment", y = "Ratio of Major to Minor Axis") +
  scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA Exposed"), values = c("skyblue2", "red2")) +
  geom_errorbar(aes(ymin= meanrat$MajMinRat-serat$MajMinRat, ymax=meanrat$MajMinRat+serat$MajMinRat),width=0.2, position=position_dodge(.9))+
  theme_classic()
## Warning: Use of `meanrat$MajMinRat` is discouraged. Use `MajMinRat` instead.

## Warning: Use of `meanrat$MajMinRat` is discouraged. Use `MajMinRat` instead.

Look at the ratio of perimeter to diameter. For a circle, this would be equal to pi

PerDiRat<- lmer(PerimDiamRat~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) +(1|CrossID), data=LarvaeDat)
## boundary (singular) fit: see ?isSingular
PerDiStep<- step(PerDiRat)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -1.2e+01
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(PerDiStep)
## Backward reduced random-effect table:
## 
##                   Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                         11 4346.6 -8671.2                         
## (1 | MaleID)               1   10 4346.6 -8673.2  0.000  1          1    
## (1 | AccTankID)            2    9 4346.6 -8675.2  0.000  1          1    
## (1 | JarSeatable)          3    8 4346.6 -8677.2  0.000  1          1    
## (1 | CrossID)              4    7 4346.6 -8679.2  0.000  1          1    
## (1 | FemaleID)             0    6 4337.6 -8663.2 17.974  1  2.239e-05 ***
## (1 | JarID)                0    6 4333.8 -8655.5 25.700  1  3.989e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated   Sum Sq  Mean Sq NumDF DenDF F value
## WCa_pHDIC:JarOmegaCalcite          0 0.011272 0.011272     1   253  5.1505
##                            Pr(>F)  
## WCa_pHDIC:JarOmegaCalcite 0.02408 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## PerimDiamRat ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) + 
##     (1 | JarID) + WCa_pHDIC:JarOmegaCalcite
#final model chosen PerimDiamRat ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) + (1 | JarID)
PerDiRatFin<- lmer(PerimDiamRat ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) + (1 | JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(PerDiRatFin) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(PerDiRatFin))
qqline(resid(PerDiRatFin)) #this looks fine

acf(resid(PerDiRatFin))

summary(PerDiRatFin) #significant effect of parent treatment, jar, and interaction, not an additive effect though. 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PerimDiamRat ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +  
##     (1 | JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: -8693.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1373 -0.5713  0.0461  0.6356  3.1915 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  JarID    (Intercept) 1.231e-04 0.011094
##  FemaleID (Intercept) 5.609e-05 0.007489
##  Residual             2.188e-03 0.046781
## Number of obs: 2698, groups:  JarID, 270; FemaleID, 15
## 
## Fixed effects:
##                             Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)                 3.084977   0.006320  33.245328 488.133  < 2e-16 ***
## WCa_pHDIC                   0.018379   0.006896  33.063404   2.665   0.0118 *  
## JarOmegaCalcite             0.049719   0.005889 253.121060   8.442 2.44e-15 ***
## WCa_pHDIC:JarOmegaCalcite  -0.014609   0.006437 252.997219  -2.269   0.0241 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC   -0.893                
## JarOmegClct -0.617  0.551         
## WC_HDIC:JOC  0.550 -0.616   -0.893
#larvae in control treatments from both parents are more circular. Larvae are shorter and squatter from OA treatments. If they were elongated, the value would be greater than pi. 

Look at the ratio of perimeter to diameter with Parent as a covariate and JarTrt as a fixed effect. For a circle, this would be equal to pi

PerDiRatJF<- lmer(PerimDiamRat~WCa_pHDIC*JarTrt+(1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) +(1|CrossID), data=LarvaeDat)
## boundary (singular) fit: see ?isSingular
PerDiStepJF<- step(PerDiRatJF)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(PerDiStepJF)
## Backward reduced random-effect table:
## 
##                   Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                         11 4346.7 -8671.4                         
## (1 | AccTankID)            1   10 4346.7 -8673.4  0.000  1          1    
## (1 | MaleID)               2    9 4346.7 -8675.4  0.000  1          1    
## (1 | JarSeatable)          3    8 4346.7 -8677.4  0.000  1          1    
## (1 | CrossID)              4    7 4346.7 -8679.4  0.000  1          1    
## (1 | FemaleID)             0    6 4337.6 -8663.3 18.124  1   2.07e-05 ***
## (1 | JarID)                0    6 4334.0 -8656.1 25.295  1   4.92e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated   Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)  
## WCa_pHDIC:JarTrt          0 0.011719 0.011719     1 252.98  5.3549 0.02147 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## PerimDiamRat ~ WCa_pHDIC + JarTrt + (1 | FemaleID) + (1 | JarID) + 
##     WCa_pHDIC:JarTrt
#final model chosen PerimDiamRat ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | JarID)
PerDiRatFinJF<- lmer(PerimDiamRat ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(PerDiRatFinJF) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(PerDiRatFinJF))
qqline(resid(PerDiRatFinJF)) #this looks fine

acf(resid(PerDiRatFinJF))

summary(PerDiRatFinJF)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PerimDiamRat ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: -8693.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1399 -0.5686  0.0462  0.6360  3.1972 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  JarID    (Intercept) 1.219e-04 0.01104 
##  FemaleID (Intercept) 5.625e-05 0.00750 
##  Residual             2.188e-03 0.04678 
## Number of obs: 2698, groups:  JarID, 270; FemaleID, 15
## 
## Fixed effects:
##                           Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)               3.139138   0.005569  20.307544 563.724  < 2e-16 ***
## WCa_pHDIC                 0.002370   0.006091  20.396908   0.389   0.7012    
## JarTrtExposed            -0.042459   0.004997 253.096458  -8.498 1.68e-15 ***
## WCa_pHDIC:JarTrtExposed   0.012634   0.005460 252.978240   2.314   0.0215 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrTrtE
## WCa_pHDIC   -0.893                
## JarTrtExpsd -0.449  0.401         
## WC_HDIC:JTE  0.402 -0.451   -0.893

Look at the ratio of perimeter to diameter with Parent and JarTrt as a fixed effects

PerDiRatFE<- lmer(PerimDiamRat~ParentTrt*JarTrt+(1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) +(1|CrossID)+(1|TrtTankID), data=LarvaeDat)
## boundary (singular) fit: see ?isSingular
PerDiStepFE<- step(PerDiRatFE)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(PerDiStepFE)
## Backward reduced random-effect table:
## 
##                   Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                         12 4346.4 -8668.8                         
## (1 | MaleID)               1   11 4346.4 -8670.8  0.000  1          1    
## (1 | AccTankID)            2   10 4346.4 -8672.8  0.000  1          1    
## (1 | TrtTankID)            3    9 4346.4 -8674.8  0.000  1          1    
## (1 | CrossID)              4    8 4346.4 -8676.8  0.000  1          1    
## (1 | JarSeatable)          5    7 4346.4 -8678.8  0.000  1          1    
## (1 | FemaleID)             0    6 4336.9 -8661.7 19.087  1  1.249e-05 ***
## (1 | JarID)                0    6 4333.9 -8655.8 25.024  1  5.663e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated   Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)  
## ParentTrt:JarTrt          0 0.012977 0.012977     1 252.99  5.9295 0.01558 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## PerimDiamRat ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 | JarID) + 
##     ParentTrt:JarTrt
#final model chosen PerimDiamRat ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | JarID)
PerDiRatFinFE<- lmer(PerimDiamRat ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 | JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(PerDiRatFinFE) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(PerDiRatFinFE))
qqline(resid(PerDiRatFinFE)) #this looks fine

acf(resid(PerDiRatFinFE))

summary(PerDiRatFinFE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PerimDiamRat ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 | JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: -8692.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1404 -0.5705  0.0466  0.6357  3.2012 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  JarID    (Intercept) 1.212e-04 0.01101 
##  FemaleID (Intercept) 5.837e-05 0.00764 
##  Residual             2.188e-03 0.04678 
## Number of obs: 2698, groups:  JarID, 270; FemaleID, 15
## 
## Fixed effects:
##                               Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                   3.141620   0.003477  20.281561 903.594  < 2e-16
## ParentTrt2800                -0.001170   0.005082  20.162026  -0.230   0.8202
## JarTrtExposed                -0.027015   0.003075 252.818361  -8.785 2.43e-16
## ParentTrt2800:JarTrtExposed  -0.010961   0.004501 252.989564  -2.435   0.0156
##                                
## (Intercept)                 ***
## ParentTrt2800                  
## JarTrtExposed               ***
## ParentTrt2800:JarTrtExposed *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.684               
## JarTrtExpsd -0.448  0.307        
## PrT2800:JTE  0.306 -0.446  -0.683

Visualize Perimeter to Diameter ratio

ggplot(LarvaeDat, aes(x=PerimDiamRat, color=JarTrt))+
  geom_histogram(alpha=0.5, position="identity") +
  facet_grid(~ParentTrt)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

meanperimdiam<- aggregate(PerimDiamRat~ParentTrt*JarTrt, data=LarvaeDat, FUN=mean)
seperimdiam<- aggregate(PerimDiamRat~ParentTrt*JarTrt, data=LarvaeDat, FUN=se)
ggplot(data = meanperimdiam, aes(x = ParentTrt, y = PerimDiamRat, fill=JarTrt)) +
  geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
  labs(x = "Parent Treatment", y = "Perimeter to Diameter Ratio") +
  scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA"), values = c("skyblue2", "red2")) +
  geom_errorbar(aes(ymin= PerimDiamRat-seperimdiam$PerimDiamRat, ymax=PerimDiamRat+seperimdiam$PerimDiamRat),width=0.2, position=position_dodge(.9))+
  geom_abline(slope=0, intercept=pi,colour="orange")+
  theme_classic()

Look at eccentricity

#look at the eccentricity data
hist(LarvaeDat$LarvaeEccentricity)#data apprear to be left skewed. 

#get the mean eccentricity for the treatments to plot
meaneccen<- aggregate(LarvaeEccentricity~ParentTrt*JarTrt, data=LarvaeDat, FUN=mean)
seeccen<- aggregate(LarvaeEccentricity~ParentTrt*JarTrt, data=LarvaeDat, FUN=se)
ggplot(data = meaneccen, aes(x = ParentTrt, y = LarvaeEccentricity, fill=JarTrt)) +
  geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
  labs(x = "Parent Treatment", y = "Larvae Eccentricity") +
  scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA"), values = c("skyblue2", "red2")) +
  geom_errorbar(aes(ymin= LarvaeEccentricity-seeccen$LarvaeEccentricity, ymax=LarvaeEccentricity+seeccen$LarvaeEccentricity),width=0.2, position=position_dodge(.9))+
  theme_classic()

jarmean<-aggregate(LarvaeEccentricity~JarTrt, data=LarvaeDat, FUN=mean)
sejar<- aggregate(LarvaeEccentricity~JarTrt, data=LarvaeDat, FUN=se) 
#plot just by JarTrt
ggplot(data = jarmean, aes(x = JarTrt, y = LarvaeEccentricity)) +
  geom_bar(stat="identity", aes(fill="gray45"), position="dodge")+
  labs(x = "Jar Treatment", y = "Larvae Eccentricity") +
  scale_fill_manual(values = c("gray45")) +
  geom_errorbar(aes(ymin= LarvaeEccentricity-sejar$LarvaeEccentricity, ymax=LarvaeEccentricity+sejar$LarvaeEccentricity),width=0.2, position=position_dodge(.9))+
  theme_classic()

larveccen<-lmer(LarvaeEccentricity~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID), data=LarvaeDat)
steppedlarveccen<- step(larveccen)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00337866 (tol = 0.002, component 1)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(steppedlarveccen)
## Backward reduced random-effect table:
## 
##                   Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                         11 3928.2 -7834.3                         
## (1 | AccTankID)            1   10 3928.2 -7836.3  0.000  1     1.0000    
## (1 | MaleID)               2    9 3928.2 -7838.3  0.003  1     0.9555    
## (1 | JarSeatable)          3    8 3928.1 -7840.2  0.083  1     0.7732    
## (1 | CrossID)              4    7 3927.4 -7840.8  1.446  1     0.2292    
## (1 | FemaleID)             0    6 3917.7 -7823.5 19.345  1  1.091e-05 ***
## (1 | JarID)                0    6 3906.9 -7801.8 40.998  1  1.524e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated   Sum Sq  Mean Sq NumDF   DenDF F value
## WCa_pHDIC:JarOmegaCalcite          1 0.010813 0.010813     1 252.181  3.6597
## WCa_pHDIC                          2 0.003051 0.003051     1  13.003  1.0327
## JarOmegaCalcite                    0 0.021731 0.021731     1 253.313  7.3556
##                             Pr(>F)   
## WCa_pHDIC:JarOmegaCalcite 0.056875 . 
## WCa_pHDIC                 0.328064   
## JarOmegaCalcite           0.007144 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## LarvaeEccentricity ~ JarOmegaCalcite + (1 | FemaleID) + (1 | 
##     JarID)
#final model only had JarTrt, female ID and jarID
larveccenfin<- lmer(LarvaeEccentricity~JarOmegaCalcite+ (1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larveccenfin) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larveccenfin))
qqline(resid(larveccenfin)) #has relatively long tails, doesn't meet assumption of normality. 

acf(resid(larveccenfin))

summary(larveccenfin) #significant effect of JarTrt
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeEccentricity ~ JarOmegaCalcite + (1 | FemaleID) + (1 |  
##     JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: -7866.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.8185 -0.4674  0.0592  0.5861  3.4958 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  JarID    (Intercept) 2.256e-04 0.015020
##  FemaleID (Intercept) 8.923e-05 0.009446
##  Residual             2.954e-03 0.054354
## Number of obs: 2698, groups:  JarID, 270; FemaleID, 15
## 
## Fixed effects:
##                  Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)     4.970e-01  3.540e-03 3.476e+01 140.388  < 2e-16 ***
## JarOmegaCalcite 8.879e-03  3.274e-03 2.533e+02   2.712  0.00714 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## JarOmegClct -0.609
#try a transformation
larveccent<-lmer(LarvaeEccentricity^3~JarOmegaCalcite+ (1|FemaleID)+(1|JarID) , data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larveccent) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larveccent))
qqline(resid(larveccent)) #has relatively long tails, doesn't meet assumption of normality. transformation helped, but it's not great

summary(larveccent) #no significant differences 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeEccentricity^3 ~ JarOmegaCalcite + (1 | FemaleID) + (1 |  
##     JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: -9778.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2074 -0.5825 -0.0178  0.5619  4.4271 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  JarID    (Intercept) 1.006e-04 0.010028
##  FemaleID (Intercept) 4.613e-05 0.006792
##  Residual             1.459e-03 0.038201
## Number of obs: 2698, groups:  JarID, 270; FemaleID, 15
## 
## Fixed effects:
##                  Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)     1.311e-01  2.488e-03 3.317e+01  52.704   <2e-16 ***
## JarOmegaCalcite 1.134e-03  2.252e-03 2.534e+02   0.504    0.615    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## JarOmegClct -0.596

Analyze eccentricity with parent as a covariate and jar as fixed effect

larveccenJF<-lmer(LarvaeEccentricity~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID), data=LarvaeDat)
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -4.3e+01
steppedlarveccenJF<- step(larveccenJF)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -4.2e+01
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
## Warning: Model failed to converge with 1 negative eigenvalue: -2.5e+01
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
## Warning: Model failed to converge with 1 negative eigenvalue: -2.5e+01
## boundary (singular) fit: see ?isSingular
print(steppedlarveccenJF)
## Backward reduced random-effect table:
## 
##                   Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                         11 3927.8 -7833.7                         
## (1 | AccTankID)            1   10 3927.8 -7835.7  0.000  1     1.0000    
## (1 | MaleID)               2    9 3927.8 -7837.7  0.000  1     0.9995    
## (1 | JarSeatable)          3    8 3927.8 -7839.6  0.087  1     0.7682    
## (1 | CrossID)              4    7 3927.1 -7840.1  1.439  1     0.2304    
## (1 | FemaleID)             0    6 3917.4 -7822.8 19.342  1  1.093e-05 ***
## (1 | JarID)                0    6 3906.6 -7801.1 41.002  1  1.521e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated    Sum Sq   Mean Sq NumDF   DenDF F value   Pr(>F)
## WCa_pHDIC:JarTrt          1 0.0108336 0.0108336     1 252.166  3.6668 0.056637
## WCa_pHDIC                 2 0.0030599 0.0030599     1  13.003  1.0357 0.327385
## JarTrt                    0 0.0216891 0.0216891     1 253.301  7.3414 0.007198
##                    
## WCa_pHDIC:JarTrt . 
## WCa_pHDIC          
## JarTrt           **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## LarvaeEccentricity ~ JarTrt + (1 | FemaleID) + (1 | JarID)
#final model
larveccenfinJF<- lmer(LarvaeEccentricity~JarTrt+ (1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larveccenfinJF) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larveccenfinJF))
qqline(resid(larveccenfinJF)) #has relatively long tails, doesn't meet assumption of normality. 

acf(resid(larveccenfinJF))

summary(larveccenfinJF)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeEccentricity ~ JarTrt + (1 | FemaleID) + (1 | JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: -7865.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.8185 -0.4705  0.0593  0.5863  3.4958 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  JarID    (Intercept) 2.257e-04 0.015023
##  FemaleID (Intercept) 8.905e-05 0.009437
##  Residual             2.954e-03 0.054354
## Number of obs: 2698, groups:  JarID, 270; FemaleID, 15
## 
## Fixed effects:
##                 Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     0.506682   0.003135  21.768315 161.632   <2e-16 ***
## JarTrtExposed  -0.007533   0.002780 253.301235  -2.709   0.0072 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## JarTrtExpsd -0.447

Analyze eccentricity with parent and jar as fixed effects

larveccenFE<-lmer(LarvaeEccentricity~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|TrtTankID)+(1|AccTankID), data=LarvaeDat)
## boundary (singular) fit: see ?isSingular
steppedlarveccenFE<- step(larveccenFE)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -1.4e+00
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00459713 (tol = 0.002, component 1)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(steppedlarveccenFE)
## Backward reduced random-effect table:
## 
##                   Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                         12 3927.5 -7831.1                         
## (1 | TrtTankID)            1   11 3927.5 -7833.1  0.000  1     1.0000    
## (1 | AccTankID)            2   10 3927.5 -7835.1  0.003  1     0.9569    
## (1 | MaleID)               3    9 3927.5 -7837.1  0.000  1     1.0000    
## (1 | JarSeatable)          4    8 3927.5 -7839.0  0.088  1     0.7664    
## (1 | CrossID)              5    7 3926.8 -7839.5  1.455  1     0.2277    
## (1 | FemaleID)             0    6 3916.8 -7821.7 19.840  1  8.418e-06 ***
## (1 | JarID)                0    6 3906.4 -7800.8 40.785  1  1.699e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated   Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)  
## ParentTrt:JarTrt          0 0.011839 0.011839     1 252.18  4.0073 0.04637 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## LarvaeEccentricity ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 | 
##     JarID) + ParentTrt:JarTrt
#final model
larveccenfinFE<- lmer(LarvaeEccentricity~ParentTrt*JarTrt+ (1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larveccenfinFE) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larveccenfinFE))
qqline(resid(larveccenfinFE)) #has relatively long tails, doesn't meet assumption of normality. 

acf(resid(larveccenfinFE))

summary(larveccenfinFE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeEccentricity ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 |  
##     JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: -7853.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.7892 -0.4764  0.0642  0.5809  3.5249 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  JarID    (Intercept) 2.193e-04 0.014809
##  FemaleID (Intercept) 9.099e-05 0.009539
##  Residual             2.954e-03 0.054355
## Number of obs: 2698, groups:  JarID, 270; FemaleID, 15
## 
## Fixed effects:
##                               Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                  5.069e-01  4.316e-03  2.012e+01 117.439  < 2e-16
## ParentTrt2800               -3.682e-04  6.309e-03  2.000e+01  -0.058 0.954041
## JarTrtExposed               -1.270e-02  3.784e-03  2.520e+02  -3.356 0.000913
## ParentTrt2800:JarTrtExposed  1.109e-02  5.538e-03  2.522e+02   2.002 0.046375
##                                
## (Intercept)                 ***
## ParentTrt2800                  
## JarTrtExposed               ***
## ParentTrt2800:JarTrtExposed *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.684               
## JarTrtExpsd -0.444  0.304        
## PrT2800:JTE  0.304 -0.442  -0.683

Test Larvae perimeter

ggplot(LarvaeDat, aes(x=LarvaePerimeterum, color=JarTrt))+
  geom_histogram(alpha=0.5, position="identity") +
  facet_grid(~ParentTrt)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#test to see if there is a significant difference
larvperim<-lmer(LarvaePerimeterum~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID), data=LarvaeDat)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0059559 (tol = 0.002, component 1)
steppedlarvperim<- step(larvperim)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00511491 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0027975 (tol = 0.002, component 1)
print(steppedlarvperim)
## Backward reduced random-effect table:
## 
##                   Eliminated npar  logLik   AIC     LRT Df Pr(>Chisq)    
## <none>                         11 -9769.7 19561                          
## (1 | AccTankID)            1   10 -9769.9 19560   0.459  1  0.4980403    
## (1 | MaleID)               2    9 -9770.9 19560   1.951  1  0.1624814    
## (1 | FemaleID)             0    8 -9773.5 19563   5.260  1  0.0218204 *  
## (1 | CrossID)              0    8 -9777.3 19571  12.689  1  0.0003677 ***
## (1 | JarSeatable)          0    8 -9773.4 19563   4.965  1  0.0258594 *  
## (1 | JarID)                0    8 -9869.4 19755 197.038  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated Sum Sq Mean Sq NumDF DenDF F value
## WCa_pHDIC:JarOmegaCalcite          0 573.37  573.37     1 221.4  8.0187
##                             Pr(>F)   
## WCa_pHDIC:JarOmegaCalcite 0.005056 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## LarvaePerimeterum ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) + 
##     (1 | CrossID) + (1 | JarSeatable) + (1 | JarID) + WCa_pHDIC:JarOmegaCalcite
#final model chosen by the step function: LarvaePerimeterum~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)
larvperimfin<- lmer(LarvaePerimeterum~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvperimfin) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvperimfin))
qqline(resid(larvperimfin)) #has relatively long tails, doesn't meet assumption of normality. 

acf(resid(larvperimfin)) #this looks fine to me.

summary(larvperimfin) #should I try to transform? Main effects and interaction are significant. 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaePerimeterum ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +  
##     (1 | CrossID) + (1 | JarSeatable) + (1 | JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: 19541.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4544 -0.5443  0.0484  0.6225  5.3183 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  JarID       (Intercept) 16.212   4.026   
##  CrossID     (Intercept)  5.489   2.343   
##  FemaleID    (Intercept)  5.599   2.366   
##  JarSeatable (Intercept)  1.107   1.052   
##  Residual                71.504   8.456   
## Number of obs: 2698, groups:  
## JarID, 270; CrossID, 45; FemaleID, 15; JarSeatable, 3
## 
## Fixed effects:
##                           Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)                196.343      2.069  24.736  94.887  < 2e-16 ***
## WCa_pHDIC                   -5.511      2.160  23.829  -2.552  0.01755 *  
## JarOmegaCalcite             37.384      1.539 221.250  24.284  < 2e-16 ***
## WCa_pHDIC:JarOmegaCalcite    4.767      1.683 221.395   2.832  0.00506 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC   -0.854                
## JarOmegClct -0.493  0.460         
## WC_HDIC:JOC  0.439 -0.514   -0.893

Analyze larvae perimeter with adult as a covariate and jar as a fixed effect

larvperimJF<-lmer(LarvaePerimeterum~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID), data=LarvaeDat)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00324995 (tol = 0.002, component 1)
steppedlarvperimJF<- step(larvperimJF)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00257919 (tol = 0.002, component 1)
print(steppedlarvperimJF)
## Backward reduced random-effect table:
## 
##                   Eliminated npar  logLik   AIC     LRT Df Pr(>Chisq)    
## <none>                         11 -9767.5 19557                          
## (1 | AccTankID)            1   10 -9767.7 19555   0.454  1  0.5003185    
## (1 | CrossID)              2    9 -9769.0 19556   2.586  1  0.1077988    
## (1 | FemaleID)             0    8 -9788.4 19593  38.693  1  4.959e-10 ***
## (1 | MaleID)               0    8 -9776.4 19569  14.715  1  0.0001251 ***
## (1 | JarSeatable)          0    8 -9772.4 19561   6.684  1  0.0097280 ** 
## (1 | JarID)                0    8 -9880.6 19777 223.234  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## WCa_pHDIC:JarTrt          0  505.5   505.5     1 240.59  7.0697 0.008366 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## LarvaePerimeterum ~ WCa_pHDIC + JarTrt + (1 | FemaleID) + (1 | 
##     MaleID) + (1 | JarSeatable) + (1 | JarID) + WCa_pHDIC:JarTrt
#final model chosen by the step function: LarvaePerimeterum~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)
larvperimfinJF<- lmer(LarvaePerimeterum~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvperimfinJF) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvperimfinJF))
qqline(resid(larvperimfinJF)) #has relatively long tails, doesn't meet assumption of normality. 

acf(resid(larvperimfinJF)) #this looks fine to me.

summary(larvperimfinJF)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaePerimeterum ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 |  
##     MaleID) + (1 | JarSeatable) + (1 | JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: 19538
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3527 -0.5461  0.0479  0.6246  5.2579 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  JarID       (Intercept) 16.880   4.108   
##  FemaleID    (Intercept)  8.381   2.895   
##  MaleID      (Intercept)  3.937   1.984   
##  JarSeatable (Intercept)  1.432   1.197   
##  Residual                71.503   8.456   
## Number of obs: 2698, groups:  
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)             237.0684     2.3822  21.0475  99.516  < 2e-16 ***
## WCa_pHDIC                -0.4725     2.5431  19.3191  -0.186  0.85453    
## JarTrtExposed           -31.8749     1.3265 240.6127 -24.029  < 2e-16 ***
## WCa_pHDIC:JarTrtExposed  -3.8550     1.4499 240.5926  -2.659  0.00837 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) WC_pHDIC JrTrtE
## WCa_pHDIC   -0.849                
## JarTrtExpsd -0.279  0.255         
## WC_HDIC:JTE  0.249 -0.287   -0.893

Analyze larvae perimeter with adult and jar as a fixed effects

larvperimFE<-lmer(LarvaePerimeterum~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|TrtTankID)+(1|AccTankID), data=LarvaeDat)
## boundary (singular) fit: see ?isSingular
steppedlarvperimFE<- step(larvperimFE)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00256065 (tol = 0.002, component 1)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00531912 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00469702 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00531912 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00281445 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.004105 (tol = 0.002, component 1)
print(steppedlarvperimFE)
## Backward reduced random-effect table:
## 
##                   Eliminated npar  logLik   AIC     LRT Df Pr(>Chisq)    
## <none>                         12 -9768.4 19561                          
## (1 | TrtTankID)            1   11 -9768.4 19559   0.000  1  0.9964895    
## (1 | AccTankID)            2   10 -9768.6 19557   0.436  1  0.5092703    
## (1 | CrossID)              3    9 -9769.9 19558   2.517  1  0.1125893    
## (1 | FemaleID)             0    8 -9788.8 19594  37.821  1  7.753e-10 ***
## (1 | MaleID)               0    8 -9777.2 19570  14.579  1  0.0001344 ***
## (1 | JarSeatable)          0    8 -9773.2 19562   6.639  1  0.0099792 ** 
## (1 | JarID)                0    8 -9882.5 19781 225.181  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## ParentTrt:JarTrt          0  425.2   425.2     1 240.57  5.9465 0.01547 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## LarvaePerimeterum ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 | 
##     MaleID) + (1 | JarSeatable) + (1 | JarID) + ParentTrt:JarTrt
#final model chosen by the step function: LarvaePerimeterum~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)
larvperimfinFE<- lmer(LarvaePerimeterum~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvperimfinFE) #seems to meet assumption of linearity and homoscedasticity

qqnorm(resid(larvperimfinFE))
qqline(resid(larvperimfinFE)) #has relatively long tails, doesn't meet assumption of normality. 

acf(resid(larvperimfinFE)) #this looks fine to me.

summary(larvperimfinFE)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaePerimeterum ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 |  
##     MaleID) + (1 | JarSeatable) + (1 | JarID)
##    Data: LarvaeDat
## 
## REML criterion at convergence: 19539.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3497 -0.5441  0.0477  0.6244  5.2637 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  JarID       (Intercept) 16.991   4.122   
##  FemaleID    (Intercept)  8.259   2.874   
##  MaleID      (Intercept)  3.925   1.981   
##  JarSeatable (Intercept)  1.430   1.196   
##  Residual                71.503   8.456   
## Number of obs: 2698, groups:  
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
## 
## Fixed effects:
##                             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                 236.3700     1.6243  16.9133 145.519   <2e-16 ***
## ParentTrt2800                 0.6594     2.0937  19.2644   0.315   0.7562    
## JarTrtExposed               -36.3906     0.8198 240.5488 -44.392   <2e-16 ***
## ParentTrt2800:JarTrtExposed   2.9251     1.1995 240.5746   2.439   0.0155 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.636               
## JarTrtExpsd -0.256  0.199        
## PrT2800:JTE  0.175 -0.288  -0.683

Visualize the perimeter data

ggplot(LarvaeDat,
       aes(fill = WCa_pHDIC, y = LarvaePerimeterum, x = JarOmegaCalcite)) +
  geom_point(aes(shape= Block, colour=WCa_pHDIC)) +
  ggtitle("Larvae area by larvae treatment")+ 
  ylab("Larvae Perimeter (um)")+
  theme_classic()

meanperim<- ddply(LarvaeDat, .(JarTrt, ParentTrt), numcolwise(mean, na.rm=T))
seperim<- ddply(LarvaeDat, .(JarTrt, ParentTrt), numcolwise(se, na.rm=T))
ggplot(data = meanperim, aes(x = ParentTrt, y = LarvaePerimeterum, fill=JarTrt)) +
  geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
  labs(x = "Parent Treatment", y = "Larvae Perimeter (um)") +
  scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA Exposed"), values = c("skyblue2", "red2")) +
  geom_errorbar(aes(ymin= LarvaePerimeterum-seperim$LarvaePerimeterum, ymax=LarvaePerimeterum+seperim$LarvaePerimeterum),width=0.2, position=position_dodge(.9))+
  theme_classic()

Analyze survival data

surv1<- lmer(RatSurv~WCa_pHDIC+ (1|FemaleID)+(1|MaleID)+(1|AccTankID), data=SurvWideDat)
## boundary (singular) fit: see ?isSingular
survstep<- step(surv1)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(survstep)
## Backward reduced random-effect table:
## 
##                 Eliminated npar logLik     AIC     LRT Df Pr(>Chisq)    
## <none>                        6 24.955 -37.911                          
## (1 | MaleID)             1    5 24.955 -39.911  0.0000  1  1.0000000    
## (1 | AccTankID)          2    4 24.312 -40.624  1.2868  1  0.2566362    
## (1 | FemaleID)           0    3 18.506 -31.012 11.6120  1  0.0006553 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated     Sum Sq    Mean Sq NumDF DenDF F value Pr(>F)
## WCa_pHDIC          1 0.00030233 0.00030233     1    13  0.0289 0.8675
## 
## Model found:
## RatSurv ~ (1 | FemaleID)
#final model chosen by step function: RatSurv~ (1|FemaleID)
survfin<- lmer(RatSurv~ 1+ (1|FemaleID), data=SurvWideDat)
plot(survfin) #this looks fine to me

qqnorm(resid(survfin))
qqline(resid(survfin)) #I think this is fine

acf(resid(survfin))

summary(survfin)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatSurv ~ 1 + (1 | FemaleID)
##    Data: SurvWideDat
## 
## REML criterion at convergence: -51.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.37540 -0.64186  0.04267  0.52410  1.83497 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  FemaleID (Intercept) 0.01122  0.1059  
##  Residual             0.01044  0.1022  
## Number of obs: 45, groups:  FemaleID, 15
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.86470    0.03131 14.00000   27.62  1.3e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#female ID and not parent treatment significantly affects survival
#I don't think a glmer is necessary this seems to meet the assumptions of lmer
survprop<- lmer(SurvChange~WCa_pHDIC+ (1|FemaleID)+(1|MaleID)+(1|AccTankID), data=SurvWideDat)
## boundary (singular) fit: see ?isSingular
survstepprop<- step(survprop)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(survstepprop)
## Backward reduced random-effect table:
## 
##                 Eliminated npar logLik     AIC     LRT Df Pr(>Chisq)    
## <none>                        6 24.955 -37.911                          
## (1 | MaleID)             1    5 24.955 -39.911  0.0000  1  1.0000000    
## (1 | AccTankID)          2    4 24.312 -40.624  1.2868  1  0.2566362    
## (1 | FemaleID)           0    3 18.506 -31.012 11.6120  1  0.0006553 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated     Sum Sq    Mean Sq NumDF DenDF F value Pr(>F)
## WCa_pHDIC          1 0.00030233 0.00030233     1    13  0.0289 0.8675
## 
## Model found:
## SurvChange ~ (1 | FemaleID)
#final model chosen by step function: SurvChange~ (1|FemaleID)
survfinprop<- lmer(SurvChange~ 1+ (1|FemaleID), data=SurvWideDat)
plot(survfinprop) #this looks fine to me

qqnorm(resid(survfinprop))
qqline(resid(survfinprop)) #I think this is fine

acf(resid(survfinprop))

summary(survfinprop)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SurvChange ~ 1 + (1 | FemaleID)
##    Data: SurvWideDat
## 
## REML criterion at convergence: -51.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.83497 -0.52410 -0.04267  0.64186  2.37540 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  FemaleID (Intercept) 0.01122  0.1059  
##  Residual             0.01044  0.1022  
## Number of obs: 45, groups:  FemaleID, 15
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.13530    0.03131 14.00000   4.322 0.000703 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#same result as survival proportion. Female ID is significant random effect but not parent treatment

Analyze survival data with parent as a fixed effect

surv1FE<- lmer(RatSurv~ParentTrt+ (1|FemaleID)+(1|MaleID)+(1|AccTankID)+(1|TrtTankID), data=SurvWideDat)
## boundary (singular) fit: see ?isSingular
survstepFE<- step(surv1FE)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(survstepFE)
## Backward reduced random-effect table:
## 
##                 Eliminated npar logLik     AIC     LRT Df Pr(>Chisq)    
## <none>                        7 24.769 -35.539                          
## (1 | MaleID)             1    6 24.769 -37.539  0.0000  1  1.0000000    
## (1 | TrtTankID)          2    5 24.769 -39.539  0.0000  1  0.9999960    
## (1 | AccTankID)          3    4 24.114 -40.228  1.3111  1  0.2522035    
## (1 | FemaleID)           0    3 18.300 -30.601 11.6267  1  0.0006501 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated    Sum Sq   Mean Sq NumDF DenDF F value Pr(>F)
## ParentTrt          1 0.0001736 0.0001736     1    13  0.0166 0.8994
## 
## Model found:
## RatSurv ~ (1 | FemaleID)
#final model chosen by step function: RatSurv~ (1|FemaleID), which is the same as survfin from above
survpropFE<- lmer(SurvChange~ParentTrt+ (1|FemaleID)+(1|MaleID)+(1|AccTankID)+(1|TrtTankID), data=SurvWideDat)
## boundary (singular) fit: see ?isSingular
survsteppropFE<- step(survpropFE)
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
print(survsteppropFE)
## Backward reduced random-effect table:
## 
##                 Eliminated npar logLik     AIC     LRT Df Pr(>Chisq)    
## <none>                        7 24.769 -35.539                          
## (1 | MaleID)             1    6 24.769 -37.539  0.0000  1  1.0000000    
## (1 | TrtTankID)          2    5 24.769 -39.539  0.0000  1  0.9999960    
## (1 | AccTankID)          3    4 24.114 -40.228  1.3111  1  0.2522035    
## (1 | FemaleID)           0    3 18.300 -30.601 11.6267  1  0.0006501 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated    Sum Sq   Mean Sq NumDF DenDF F value Pr(>F)
## ParentTrt          1 0.0001736 0.0001736     1    13  0.0166 0.8994
## 
## Model found:
## SurvChange ~ (1 | FemaleID)
#final model chosen by step function: SurvChange~ (1|FemaleID), which is the same as survfinprop

Visualize larvae survival

par(mfcol=c(1,1))
boxplot(RatSurv~as.factor(ParentTrt), data=SurvWideDat)

boxplot(SurvChange~as.factor(ParentTrt), data= SurvWideDat)

#make reaction norms for survival color will be parent treatment. X will be jar treatment. Each cross will be its own line
ggplot(SurvDat,aes(x=JarTrt, y=LarvaeSurvived, color = ParentTrt, group= CrossID)) +
    geom_point(aes(color=ParentTrt))+
  geom_line(aes(color=ParentTrt, group=CrossID))+
  geom_errorbar(aes(ymin= LarvaeSurvived-SELarvaeSurvived, ymax=LarvaeSurvived+SELarvaeSurvived),width=0.05, position=position_dodge(0))+
    labs(x = "Larvae Treatment", y = "Number of Larvae Survived") +
  scale_color_manual(values = c("skyblue2", "red2")) +
  theme_classic()

ggplot(SurvWideDat,
       aes(y = SurvChange, x = ParentTrt)) +
  geom_point() +
  ggtitle("Larvae survival by parent treatment")+ 
  ylab("Proportion change in the number of larvae that survived")+
  theme_classic()

meansurv<- ddply(SurvWideDat, .(ParentTrt), numcolwise(mean, na.rm=T))
sesurv<- ddply(SurvWideDat, .(ParentTrt), numcolwise(se, na.rm=T))
#make a shorter name for the CrossIDs in SurvWideDat
SurvWideDat$Cross<- substr(SurvWideDat$CrossID, 1, 9)
ggplot(data = SurvWideDat, aes(x = Cross, y = RatSurv, fill=ParentTrt)) +
  geom_bar(stat="identity", position="dodge")+
  labs(x = "Cross ID", y = "Proportion of Larvae that Survived") +
  scale_fill_manual(name = "Parent Treatment", labels = c("Control", "OA Exposed"), values = c("skyblue2", "red2")) +
  theme(axis.text.x=element_text(angle=90))

ggplot(data = meansurv, aes(x = ParentTrt, y = RatSurv)) +
  geom_bar(stat="identity", position="dodge")+
  labs(x = "Parent Treatment", y = "Proportion of Larvae that Survived")+
  geom_errorbar(aes(ymin= RatSurv-sesurv$RatSurv, ymax=RatSurv+sesurv$RatSurv),width=0.2, position=position_dodge(.9))+
  theme_classic()

Cilia extrusion statistics

#Plot cilia extrusion by size to see if there is a relationship (see Waldbusser et al 2015)
#to do this Elise looked at a subset of Larvae for each jar treatment. 
par(mfcol=c(1,1))
testcil<- read.csv("../data/LarvMorphTest.csv")
testcil$CilExtruded<- as.factor(testcil$CilExtruded)
testcilsub<- subset(testcil, CilExtruded !="")
testcilsub$Flag<- as.factor(testcilsub$Flag)
testcilsub<- subset(testcilsub, Flag !="TRUE")
#take a look at the data
ggplot(testcilsub,
       aes(fill = JarTrt, y = MaxFeretDiamum, x = CilExtruded)) +
  geom_point(aes( shape=JarTrt, colour=JarTrt)) +
  ggtitle("Larvae cilia extrusion by size")+ 
  ylab("Larvae Diameter (um)")+
  facet_grid(~JarTrt)+
  theme_classic()

#test for a relationship between size and cilia extruded within each jar treatment type. 
exp<- subset(testcilsub, JarTrt=="Exposed")
con<- subset(testcilsub, JarTrt=="Control")

expcil1<-glm(CilExtruded~MaxFeretDiamum, data=exp, family=binomial)
summary(expcil1) #no significant relationship between size and cilia extrusion
## 
## Call:
## glm(formula = CilExtruded ~ MaxFeretDiamum, family = binomial, 
##     data = exp)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7172   0.7201   0.7444   0.7527   0.7725  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)
## (Intercept)    0.566735   3.046585   0.186    0.852
## MaxFeretDiamum 0.008725   0.047015   0.186    0.853
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 182.22  on 163  degrees of freedom
## Residual deviance: 182.18  on 162  degrees of freedom
## AIC: 186.18
## 
## Number of Fisher Scoring iterations: 4
boxplot(MaxFeretDiamum~CilExtruded,  data=exp)

#now test the control jars for relationship between size and cilia extrusion. 
concil1<-glm(CilExtruded~MaxFeretDiamum, data=con, family=binomial)
summary(concil1) #looks like there is a significant relationship between size and cilia extrusion
## 
## Call:
## glm(formula = CilExtruded ~ MaxFeretDiamum, family = binomial, 
##     data = con)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.9324  -0.3705  -0.3098  -0.2538   2.8520  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)  
## (Intercept)     14.8599     7.5190   1.976   0.0481 *
## MaxFeretDiamum  -0.2376     0.1026  -2.316   0.0206 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 82.805  on 179  degrees of freedom
## Residual deviance: 77.639  on 178  degrees of freedom
## AIC: 81.639
## 
## Number of Fisher Scoring iterations: 6
boxplot(MaxFeretDiamum~CilExtruded,  data=con)

Export Figures

plots.dir.path<- list.files(tempdir(), pattern="rs-graphics", full.names= TRUE)
plots.png.paths<- list.files(plots.dir.path, pattern=".png", full.names=TRUE)
file.copy(from=plots.png.paths, to="../results")
## logical(0)

Look at weight per larvae

#we can check these data in a few ways: 1) look to see if replicate jars have similar weights; 2) check to see if first and second filters measured come up with similar values for weight per larva
#Start with 1) 
#start by using all data: 
boxplot(MeanWtPerLarvae~CrossID*JarTrt, data=LarvByJarDat)

#looks like some of the crosses have high variance. could we just eliminate those? 
meanwtall<- ddply(LarvByJarDat, .(CrossID,ParentTrt, JarTrt, FemaleID, MaleID), numcolwise(mean, na.rm=T))
sewtall<- ddply(LarvByJarDat, .(CrossID, ParentTrt, JarTrt, FemaleID, MaleID), numcolwise(se, na.rm=T))
ggplot(data = meanwtall, aes(x = CrossID, y = MeanWtPerLarvae)) + 
  geom_errorbar(aes(ymin= MeanWtPerLarvae-sewtall$MeanWtPerLarvae, ymax=MeanWtPerLarvae+sewtall$MeanWtPerLarvae),width=0.2)+
  geom_point(aes(colour=JarTrt)) + 
  labs(x = "Cross ID", y = "Mean Weight Per Larvae ug") + 
   theme_classic()

#let's see if removing filters with fewer than 100 individuals helps this out and those with crystals.
subdat<- LarvByJarDat
subdat$MeanWtPerLarvae<- ifelse(subdat$F1Crystals=="FALSE", paste(subdat$MeanWtPerLarvae), paste("NA"))
subdat$MeanWtPerLarvae<- ifelse(subdat$F1LarvaeCount>100, paste(subdat$MeanWtPerLarvae), paste("NA"))
subdat$MeanWtPerLarvae<- ifelse(subdat$JarID!="J004_B1", paste(subdat$MeanWtPerLarvae), paste("NA"))#this is the only F2 filter after subsetting that has <100 larvae
subdat$JarSeatable<- as.factor(subdat$JarSeatable)
subdat$MeanWtPerLarvae<- as.numeric(subdat$MeanWtPerLarvae)
## Warning: NAs introduced by coercion
subdat$JarTrt<- as.factor(subdat$JarTrt)
subdat$LarvDensity<- subdat$MeanWtPerLarvae/subdat$LarvaeAreaum2
boxplot(MeanWtPerLarvae~CrossID*JarTrt, data=subdat)

#looks like some of the crosses have high variance. could we just eliminate those? 
submeanwtall<- ddply(subdat, .(Block, ParentTrt, JarTrt, CrossID, FemaleID, MaleID), numcolwise(mean, na.rm=T))
subsewtall<- ddply(subdat, .(Block, ParentTrt, JarTrt, CrossID, FemaleID, MaleID), numcolwise(se, na.rm=T))
ggplot(data = submeanwtall, aes(x = CrossID, y = MeanWtPerLarvae)) + 
  geom_errorbar(aes(ymin= MeanWtPerLarvae-subsewtall$MeanWtPerLarvae, ymax=MeanWtPerLarvae+subsewtall$MeanWtPerLarvae),width=0.2)+
  geom_point(aes(colour=JarTrt)) + 
  labs(x = "Cross ID", y = "Mean Weight Per Larvae ug") + 
  theme(axis.text.x = element_text(angle = 90))
## Warning: Removed 4 rows containing missing values (geom_point).

#there are some crosses that only have a single point. Those should be removed. Do this by looking at se
#first join SE to the dataset
colnames(subsewtall)<- paste("SE", colnames(subsewtall), sep="")
#merge the standard error data to the SurvDat dataframe
SubWt<- merge(submeanwtall, subsewtall, by.x= c("CrossID", "JarTrt"), by.y= c("SECrossID", "SEJarTrt"))
relcrosses<- SubWt
relcrosses$MeanWtPerLarvae<- ifelse(relcrosses$SEMeanWtPerLarvae>0, paste(relcrosses$MeanWtPerLarvae), paste("NA"))
relcrosses$MeanWtPerLarvae<- as.numeric(relcrosses$MeanWtPerLarvae)
ggplot(data = relcrosses, aes(x = CrossID, y = MeanWtPerLarvae)) + 
  geom_errorbar(aes(ymin= MeanWtPerLarvae-SEMeanWtPerLarvae, ymax=MeanWtPerLarvae+SEMeanWtPerLarvae),width=0.2)+
  geom_point(aes(colour=JarTrt)) + 
  labs(x = "Cross ID", y = "Mean Weight Per Larvae ug") + 
  theme(axis.text.x = element_text(angle = 90))
## Warning: Removed 22 rows containing missing values (geom_point).

boxplot(SEMeanWtPerLarvae~JarTrt, data= relcrosses)

#I don't know how to choose which ones to include. For now, let's say if the SE is greater than mean(MeanWtPerLarvae) then I will remove it
MeanWtVal<- mean(relcrosses$MeanWtPerLarvae, na.rm=TRUE)
relcrosses$MeanWtPerLarvae<- ifelse(relcrosses$SEMeanWtPerLarvae<= MeanWtVal, paste(relcrosses$MeanWtPerLarvae), paste("NA"))
relcrosses$MeanWtPerLarvae<- as.numeric(relcrosses$MeanWtPerLarvae)
## Warning: NAs introduced by coercion
boxplot(MeanWtPerLarvae~ParentTrt*as.factor(JarTrt), data=relcrosses)

ggplot(data = relcrosses, aes(x = CrossID, y = LarvDensity)) + 
  geom_errorbar(aes(ymin= LarvDensity-SELarvDensity, ymax=LarvDensity+SELarvDensity),width=0.2)+
  geom_point(aes(colour=JarTrt)) + 
  labs(x = "Cross ID", y = "Mean Density Per Larvae ug/um2") + 
  theme(axis.text.x = element_text(angle = 90))
## Warning: Removed 4 rows containing missing values (geom_point).

#now try to create a linear model 
Density1<- lm(LarvDensity~ParentTrt*as.factor(JarTrt), data=relcrosses)
summary(Density1)
## 
## Call:
## lm(formula = LarvDensity ~ ParentTrt * as.factor(JarTrt), data = relcrosses)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.607e-10 -2.152e-10 -5.552e-11  1.227e-10  1.144e-09 
## 
## Coefficients:
##                                          Estimate Std. Error t value Pr(>|t|)
## (Intercept)                             4.929e-10  7.358e-11   6.698 2.45e-09
## ParentTrt2800                          -1.336e-11  1.091e-10  -0.122   0.9029
## as.factor(JarTrt)Exposed                2.195e-10  1.064e-10   2.063   0.0423
## ParentTrt2800:as.factor(JarTrt)Exposed  1.515e-10  1.559e-10   0.972   0.3340
##                                           
## (Intercept)                            ***
## ParentTrt2800                             
## as.factor(JarTrt)Exposed               *  
## ParentTrt2800:as.factor(JarTrt)Exposed    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.605e-10 on 82 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.1598, Adjusted R-squared:  0.129 
## F-statistic: 5.198 on 3 and 82 DF,  p-value: 0.002452
#but from the plot perhaps something is relevant for the difference going from control to exposed conditions. 
#use spread function
relcrossestowide<- unite(relcrosses, Value, c(MeanWtPerLarvae, SEMeanWtPerLarvae, LarvDensity, SELarvDensity), sep="_")
relcrossestowide<- subset(relcrossestowide, select=c("CrossID", "Value", "JarTrt", "ParentTrt", "MaleID", "FemaleID"))
#use spread function in tidyr to make the column in long form vs. wide form
relwide<- spread(relcrossestowide, JarTrt, value= Value)
#now split the mean and SE apart again
relwide<- separate(relwide, Control, into= c("MeanWtPerLarvaeCon", "SEMeanWtPerLarvaeCon", "LarvDensityCon", "SELarvDensityCon"), sep="_")
relwide<- separate(relwide, Exposed, into= c("MeanWtPerLarvaeExp", "SEMeanWtPerLarvaeExp", "LarvDensityExp", "SELarvDensityExp"), sep="_")
#make the means and ses numeric
relwide$MeanWtPerLarvaeCon<- as.numeric(relwide$MeanWtPerLarvaeCon)
## Warning: NAs introduced by coercion
relwide$MeanWtPerLarvaeExp<- as.numeric(relwide$MeanWtPerLarvaeExp)
## Warning: NAs introduced by coercion
relwide$SEMeanWtPerLarvaeCon<- as.numeric(relwide$SEMeanWtPerLarvaeCon)
## Warning: NAs introduced by coercion
relwide$SEMeanWtPerLarvaeExp<- as.numeric(relwide$SEMeanWtPerLarvaeExp)
## Warning: NAs introduced by coercion
relwide$LarvDensityCon<- as.numeric(relwide$LarvDensityCon)
relwide$SELarvDensityCon<- as.numeric(relwide$SELarvDensityCon)
## Warning: NAs introduced by coercion
relwide$LarvDensityExp<- as.numeric(relwide$LarvDensityExp)
relwide$SELarvDensityExp<- as.numeric(relwide$SELarvDensityExp)
## Warning: NAs introduced by coercion
#get the difference from control to OA
relwide$WtDiff<- (relwide$MeanWtPerLarvaeExp-relwide$MeanWtPerLarvaeCon)
plot(WtDiff~ParentTrt, data=relwide)

boxplot(WtDiff~FemaleID, data=relwide)
abline(a=0, b=0)

relwide$DensDiff<- (relwide$LarvDensityExp-relwide$LarvDensityCon)
plot(DensDiff~ParentTrt, data=relwide)

#now try a linear model. 
Dens1<- lm(DensDiff~ParentTrt, data=relwide)
summary(Dens1)
## 
## Call:
## lm(formula = DensDiff ~ ParentTrt, data = relwide)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -5.901e-10 -1.801e-10 -1.457e-11  1.668e-10  7.025e-10 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   2.007e-10  6.006e-11   3.342  0.00185 **
## ParentTrt2800 1.184e-10  8.822e-11   1.342  0.18751   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.817e-10 on 39 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.04411,    Adjusted R-squared:  0.0196 
## F-statistic:   1.8 on 1 and 39 DF,  p-value: 0.1875
#Look at the difference between F1 and F2 measurements, so subset the data for only those with F2 measurements
#get the mean of weight per larva 
mean(LarvByJarDat$MeanWtPerLarvae)
## [1] NA
F1F2<- subset(LarvByJarDat, F2WtPerLarvae >0)
F1F2<- subset(F1F2, F1WtPerLarvae>0)
F1F2$JarTrt<- as.factor(F1F2$JarTrt)
F1F2$F1FilterID<- as.factor(F1F2$F1FilterID)
F1F2$F2FilterID<- as.factor(F1F2$F2FilterID)
F1F2$ParentTrt<- as.factor(F1F2$ParentTrt)
F1F2$CrossID<- as.factor(F1F2$CrossID)
F1F2$TrtTankID<- as.factor(F1F2$TrtTankID)
F1F2$JarSeatable<- as.factor(F1F2$JarSeatable)
F1F2$MaleID<- as.factor(F1F2$MaleID)
F1F2$FemaleID<- as.factor(F1F2$FemaleID)
F1F2$JarID<-droplevels(F1F2)$JarID
F1F2<- subset(F1F2, select=c("JarID", "F1WtPerLarvae", "F1LarvaeCount", "F2WtPerLarvae", "F2LarvaeCount", "F1FilterID", "F2FilterID", "JarTrt", "ParentTrt", "CrossID", "TrtTankID", "JarSeatable", "MaleID", "FemaleID", "F1Crystals", "MeanWtPerLarvae"))
F1F2$LarvWtDiff<- abs(F1F2$F1WtPerLarvae-F1F2$F2WtPerLarvae)
F1F2$LarvWtDiffProp<- (F1F2$LarvWtDiff)/ (F1F2$F2WtPerLarvae)
F1F2$TotalLarvCount<- F1F2$F1LarvaeCount+F1F2$F2LarvaeCount
F1F2$LarvaeCountDiff<- abs(F1F2$F1LarvaeCount- F1F2$F2WtPerLarvae)
F1F2sub<- subset(F1F2, F1LarvaeCount>100)
F1F2sub<- subset(F1F2sub, F2LarvaeCount>100)
F1F2sub<- subset(F1F2sub, F1Crystals=="FALSE")
F1F2sub$JarID<-droplevels(F1F2sub)$JarID
#Now look at the differences.In weights
mean(F1F2sub$F1WtPerLarvae)
## [1] 2.213654e-06
mean(F1F2sub$F2WtPerLarvae)
## [1] 1.594731e-06
plot(LarvWtDiff~TotalLarvCount, data=F1F2sub, col="black", pch=19) 
abline(a= 1.59e-06, b=0, col="blue")
abline(a=(2* 1.59e-06), b=0, col="red")

plot(F1LarvaeCount~LarvWtDiffProp, data=F1F2sub, pch=19)
points(F2LarvaeCount~LarvWtDiffProp, data=F1F2sub, pch=19, col="orange")

#maybe I should just remove those with larvae counts <200
F200sub<- subset(F1F2, F1LarvaeCount>200)
F200sub<- subset(F200sub, F2LarvaeCount>200)
F200sub$JarID<-droplevels(F200sub)$JarID
mean(F200sub$F1WtPerLarvae)
## [1] 1.91e-06
mean(F200sub$F2WtPerLarvae)
## [1] 1.256176e-06
plot(LarvWtDiff~TotalLarvCount, data=F200sub, col="black", pch=19) 

plot(F1LarvaeCount~LarvWtDiffProp, data=F200sub, pch=19)
points(F2LarvaeCount~LarvWtDiffProp, data=F200sub, pch=19, col="orange")

plot(F1WtPerLarvae~JarID, data=F200sub)
points(F2WtPerLarvae~JarID, data=F200sub)
points(LarvWtDiff~JarID, data=F200sub, col="red", pch=19) 

boxplot(LarvWtDiffProp~ParentTrt*JarTrt, data=F1F2sub)

#I'm going to just look at the jars that had proportion difference <1
Less1<- subset(F1F2, LarvWtDiffProp<1)
boxplot(MeanWtPerLarvae~ParentTrt*JarTrt, data=Less1)

boxplot(LarvWtDiffProp~ParentTrt*JarTrt, data=Less1)

Wtdiff<- lm(LarvWtDiff~TotalLarvCount, data=F1F2)
par(mfcol=c(2,2))
plot(Wtdiff)

par(mfcol=c(1,1))
summary(Wtdiff)
## 
## Call:
## lm(formula = LarvWtDiff ~ TotalLarvCount, data = F1F2)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.642e-06 -4.206e-06 -2.710e-06  1.996e-06  3.599e-05 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     9.669e-06  2.454e-06    3.94 0.000347 ***
## TotalLarvCount -6.907e-09  3.153e-09   -2.19 0.034877 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.982e-06 on 37 degrees of freedom
## Multiple R-squared:  0.1148, Adjusted R-squared:  0.09086 
## F-statistic: 4.798 on 1 and 37 DF,  p-value: 0.03488
#Use the Wt diff model to decide which filters to include
#first get the average larva weight overall: 
mean(F1F2$MeanWtPerLarvae)
## [1] 3.874513e-06
boxplot(F1WtPerLarvae~ParentTrt+JarTrt, data=Calc)
## Error in eval(m$data, parent.frame()): object 'Calc' not found
#I don't know whether or not to trust these data. Let's look at F1 and F2 together
Calcmean<- ddply(Calc, .(ParentTrt, JarTrt), numcolwise(mean, na.rm=TRUE))
## Error in empty(.data): object 'Calc' not found
Calcse<- ddply(Calc, .(ParentTrt, JarTrt), numcolwise(se, na.rm=TRUE))
## Error in empty(.data): object 'Calc' not found
ggplot(data = Calcmean, aes(x = ParentTrt , y = F1WtPerLarvae, fill=JarTrt)) +
  geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
  labs(x = "Parent Treatment", y = "Weight per Larvae") +
  scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA Exposed"), values = c("skyblue2", "red2")) +
  geom_errorbar(aes(ymin= F1WtPerLarvae-Calcse$F1WtPerLarvae, ymax=F1WtPerLarvae+Calcse$F1WtPerLarvae),width=0.2, position=position_dodge(.9))+
  theme_classic()
## Error in ggplot(data = Calcmean, aes(x = ParentTrt, y = F1WtPerLarvae, : object 'Calcmean' not found
calcexp<- subset(Calc, JarTrt=="Exposed")
## Error in subset(Calc, JarTrt == "Exposed"): object 'Calc' not found
calcexp$JarTrt<- as.factor(calcexp$JarTrt)
## Error in is.factor(x): object 'calcexp' not found
calccon<- subset(Calc, JarTrt=="Control")
## Error in subset(Calc, JarTrt == "Control"): object 'Calc' not found
calccon$JarTrt<- as.factor(calccon$JarTrt)
## Error in is.factor(x): object 'calccon' not found
plot(F1WtPerLarvae~F1LarvaeCount, data=calccon, col="blue", pch=19)
## Error in eval(m$data, eframe): object 'calccon' not found
points(F1WtPerLarvae~F1LarvaeCount, data=calcexp, col="red", pch=19)
## Error in eval(m$data, eframe): object 'calcexp' not found
summary(lm(F1WtPerLarvae~F1LarvaeCount+JarTrt+ParentTrt, data=Calc))
## Error in is.data.frame(data): object 'Calc' not found